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Abstract 


The  formulation  of  stable  time-stepping  algorithms  for  dynamic  contact  problems,  both 
frictionless  and  frictional,  is  presented.  Special  attention  is  given  to  the  properties  of  the 
underlying  continuum  problem  to  serve  as  guidelines  for  the  development  of  the  algorithms. 
The  proposed  method  conserves  linear  and  angular  momenta,  and,  in  the  frictionless  case, 
conserves  the  energy  by  means  of  a  restoration  potential.  Coulomb’s  friction  law  is  used  to 
model  the  friction  phenomenon;  the  scheme  presented  herein  is  unconditionally  dissipative, 
just  as  the  physical  system  is. 

The  scheme  has  been  enhanced  by  the  enforcement  of  a  constraint  on  the  velocities,  in 
addition  to  the  unilateral  (impenetrability)  constraint  imposed  on  the  displacements;  this 
enhancement  does  not  disturb  the  conservation/restoration  properties.  Numerical  dissipation 
may  also  be  added  to  stabilize  the  scheme  for  problems  with  high  frequency  energy  modes. 

A  multibody  implementation  is  presented  to  show  the  versatility  of  the  algorithm.  In  this 
implementation,  the  contact  detection  scheme  includes  an  efficient  sorting  procedure  which 
makes  large  scale  simulations  possible. 

Lastly,  various  numerical  examples  show  the  stability  and  robustness  of  the  scheme. 
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Chapter  1 


Introduction 


1.1  Motivation 

Contact  or  impact  situations  are  present  in  many  engineering  applications,  such  as  metal 
forming  processes  and  crash- worthiness  testing.  The  modeling  of  contact  problems  is  a  highly 
nonlinear  situation  due  to  the  fact  that  it  is  a  unilaterally  constrained  problem.  The  intro¬ 
duction  of  friction  between  the  contacting  bodies  is  another  important  source  of  nonlinearity. 
The  frictionless  case  is  in  fact  a  unilaterally  constrained  Hamiltonian  system  giving  rise  to 
many  conserving  properties.  The  challenge  is  for  the  algorithmic  schemes  to  simulate  these 
conservation  properties  accurately. 

To  circumvent  the  difficulties  associated  with  multiple  nonlinearities,  many  explicit  schemes 
have  been  developed  in  the  past.  But  one  of  the  main  drawbacks  of  this  type  of  scheme 
becomes  evident  when  modeling  dynamic  problems.  Explicit  schemes  have  limited  stability 
properties,  sometimes  leading  to  poor  enforcement  of  the  constraints.  On  the  other  hand, 
some  implicit  schemes  are  stable  in  a  linear  regime,  but  may  loose  their  stability  properties 
in  a  nonlinear  problem,  giving  rise  to  a  non-physical  energy  increase  in  the  system.  Char¬ 
acteristic  examples  are  the  trapezoidal  and  mid-point  rules,  which  are  energy  conserving 
schemes  in  the  linear  regime  but  show  significant  increase  in  the  energy  of  the  system  for 
nonlinear  problems,  and  may  even  result  in  numerical  blow-ups. 
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1.2.  BACKGROUND 


In  addition  to  the  impenetrability  constraint,  the  system  must  be  constrained  so  that  the  rate 
of  separation  during  contact  is  zero.  In  fact,  higher  order  time  derivatives  of  this  measure 
of  separation  are  also  required  to  vanish.  A  poor  enforcement  of  these  constraints  yields 
oscillations  between  contact  and  release  states  which  damage  the  robustness  of  a  contact 
scheme. 

The  source  of  the  oscillatory  behaviour  is  the  inability  of  a  space  and  time  discretized  for¬ 
mulation  to  represent  the  shock  wave  that  reverses  the  velocity  of  the  contact  point  at  the 
instant  of  contact.  The  ineffectiveness  of  representing  these  high  frequency  energy  modes, 
known  also  as  the  Gibbs  phenomenon,  is  a  well  knoAvn  source  of  difficulty  in  impact  algo¬ 
rithms  where  the  short  term  behaviour  is  of  interest. 

The  simulation  of  contact  problems  also  involves  contact  detection  algorithms  which  may 
affect  the  overall  efficiency  of  the  scheme,  especially  in  cases  of  many  body  problems.  A 
robust  contact  scheme  should  be  able  to  resolve  multiple  collisions  among  multiple  bodies. 
All  of  these  considerations  play  an  important  role  in  the  design  of  both  the  contact  detection 
algorithm  and  the  overall  contact  formulation. 


1.2  Background 

Since  the  development  of  the  finite  element  method  in  the  late  1950’s,  numerical  solutions  to 
contact  problems  have  been  investigated  intensively  by  many  researchers.  Earlier  work  done 
by  Francavilla  &  ZIENKIEWICZ  [10]  on  a  flexibility  approach  and  by  Hughes  et  al[15] 
on  the  use  of  Lagrange  multiplier  methods,  contributed  to  the  development  of  robust  finite 
element  methods.  A  fairly  comprehensive  overview  of  the  numerical  methodologies  used  to 
solve  the  quasi-static  contact  problem  can  be  found  in  Zhong  Mackerle[38]. 

Descriptions  of  finite  elements  methods  for  dynamic  contact  problems  may  be  found  in 
Belytschko  &  Neal[5]  who  developed  pinball  methodologies;  in  Carpenter  et  al[7], 
who  worked  on  the  development  of  Lagrange  multipliers  for  the  enforcement  of  unilateral 
constraints;  and  in  Hallquist  et  al[12],  who  developed  the  concept  of  master  and  slave 
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methodology.  All  of  these  formulations  involve  explicit  integrators. 

A  comprehensive  review  of  both  the  frictionless  and  frictional  dynamic  contact  problems  is 
presented  in  Kikuchi  &  Oden[16]. 

Enforcement  of  higher  order  constraints  (i.e. vanishing  of  rate  of  separation)  in  addition  to 
the  impenetrability  constraint  has  been  researched  by  Lee  [21]  and  Taylor  &  Papadopou- 
los[33]. 

When  dealing  with  nonlinear  dynamical  systems,  stability  issues  become  critical  in  the  sim¬ 
ulation  of  dynamic  contact  problems.  In  the  area  of  elastodynamics,  the  consideration  of 
energy /momentum  conserving  algorithms  leading  to  stable  schemes,  is  described  in  SiMO  & 
Tarnow  [30]  and  Simo  ET  al[31].  Examples  of  the  work  concentrating  on  stability  and 
conservation  properties  include  Munjiza  et  al[23]  and  Armero  &  Petocz[1]. 

In  the  area  of  frictional  contact  problems,  formulations  using  both  penalty  regularization  and 
augmented  Lagrangians  has  been  published  by  Laursen  &  Simo[19],  Oden  &  Martins 
[25]  and  Wriggers  et  al[37]. 

Contact  detection  algorithms  for  multibody  problems  have  been  developed  primarily  by 
those  in  the  discrete  element  method  research  community,  who  deal  predominantly  with  rigid 
bodies.  However,  many  ideas  can  be  extracted  from  their  work  to  be  used  in  conjunction  with 
the  finite  element  method.  A  state-of-the-art  review  is  included  in  Williams  &  O ’CONNOR 
[35].  Some  examples  of  contact  detection  algorithms  using  tree  data  structures  are  presented 
in  Munjiza  et  al[24]  and  Bonet  &  Peraire[6]. 


1.3  Goals 

The  main  goal  of  this  work  is  to  develop  stable  implicit  time-stepping  algorithms  for  dynamic 
contact  problems.  We  have  analyzed  the  physical  properties  of  the  continuum  dynamical 
system  and  developed  our  algorithm  following  these  guidelines. 

In  a  general  contact  problem  among  elastic  bodies,  in  the  absence  of  external  forces  and 
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imposed  displacements,  linear  and  angular  momenta  are  conserved.  The  proposed  scheme 
inherits  these  conservation  properties  by  construction.  Furthermore,  in  the  frictionless  case, 
the  total  energy  of  the  system  is  also  conserved.  The  algorithm  presented  in  this  work 
accomplishes  the  restoration  of  the  total  energy  after  the  release  of  the  contacting  bodies. 
That  is,  the  penalty  regularization  potential  used  to  enforce  the  unilateral  constraint  is  also 
used  to  restore  the  energy  that  was  taken  away  from  the  system  during  contact.  The  energy 
restoration  scheme  is  second  order  accurate  and  is  unconditionally  (energy)  stable  without 
relying  on  artificial  numerical  dissipation. 

During  persistent  contact,  in  which  two  bodies  may  remain  in  contact  for  long  periods  of 
time,  we  have  additional  constraints  which  force  the  rate  of  separation  between  the  contact 
surfaces  to  be  zero. 

As  mentioned  above,  the  presence  of  high  frequency  modes  in  impact  problems  tests  the 
robustness  of  many  contact  schemes.  Numerical  dissipation  is  added  to  the  developed  scheme 
in  order  to  circumvent  the  difficulties  inherent  in  this  problem.  The  amount  of  dissipation 
introduced  into  the  scheme  can  be  regulated  to  satisfy  the  requirements  of  each  particular 
problem. 

Friction  phenomena  are  purely  dissipative  physical  phenomena;  hence  our  goal  is  to  develop 
a  frictional  scheme  which  can  guarantee  positive  energy  dissipation  under  any  circumstances. 
The  proposed  scheme  is  unconditionally  dissipative. 

The  contact  algorithm  has  to  be  easy  to  implement  within  a  multibody  contact  formulation. 
The  contact  detection  scheme  is  modified  by  means  of  object-oriented  programming  and 
intelligent  database  structures,  to  accommodate  the  new  contact  detection  schemes. 

By  a  series  of  numerical  examples,  we  show  the  different  properties  of  the  contact  scheme. 
One  dimensional  benchmark  problems  and  two  dimensional  quasi  static  and  dynamical  ex¬ 
amples  clearly  show  the  superior  stability  and  conservation  properties  the  proposed  scheme 
over  conventional  schemes. 
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1.4  Outline 

In  Chapter  2  we  shall  the  describe  the  problem  of  contact  between  elastic  bodies  and  explain 
the  notation  used  in  this  work.  We  state  the  governing  equations  for  both  the  frictionless 
and  the  frictional  case.  In  Chapter  3,  we  develop  the  weak  formulation  for  the  previously 
stated  strong  form  of  the  equations  and  study  its  properties.  In  Chapter  4,  we  describe 
the  finite  element  implementation  and  and  time-stepping  procedure  to  solve  the  dynamical 
problem.  Contact  detection  in  a  space  discretized  setting  will  also  be  discussed.  We  describe 
the  the  algorithms  developed  for  the  frictionless  contact  problem  and  for  the  frictional  case 
in  Chapters  5  and  6,  respectively,  and  we  study  their  stability  properties.  In  Chapter  7,  we 
go  on  to  describe  the  multibody  implementation  of  the  contact  schemes  developed  in  this 
work.  In  Chapter  8,  we  present  various  numerical  examples  in  one  and  two  dimensions  to 
show  the  performance  of  our  schemes.  We  draw  conclusions  and  propose  some  future  work 
in  Chapter  9. 


Chapter  2 


Problem  Definition 


2.1  Introduction 

In  this  chapter,  we  develop  the  equations  that  define  the  contact  problem  and  the  nota¬ 
tion  required  to  interpret  them.  The  equations  describe  the  continuum  problem  of  contact 
between  two  or  more  bodies  with  large  deformations,  for  both  frictionless  and  frictional  cases. 

The  equations  are  defined  using  the  Lagrangian  description,  which  is  the  most  prevalent 
framework  in  computational  solid  mechanics.  For  the  sake  of  generality,  the  equations  and 
subsequent  algorithmic  solutions  will  be  developed  in  a  dynamic  formulation  but  the  work 
presented  herein  is  also  valid  for  the  quasi-static  case. 

Contact  is  a  unilaterally  constrained  problem  in  which  we  impose  that  two  bodies  may  not 
penetrate  each  other.  Relative  tangent  motion  may  be  present  and  friction  may  be  included 
in  the  statement  of  the  problem.  Many  models  have  been  proposed  to  describe  the  friction 
phenomena,  however  in  the  context  of  the  finite  element  method  mo;t  rely  on  the  classical 
law  of  Coulomb.  A  description  of  other  friction  models  may  be  found  in  Oden  &  Martins 
[25], 

We  introduce  basic  notation  in  Section  2,  describe  the  governing  equations  in  their  local 
form  in  Section  3  and  describe  the  frictional  problem  in  Section  4. 
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2.2  Notation 

We  describe  the  contact  problem  between  two  elastic  bodies;  special  cases  such  as  self-contact 
are  excluded  from  the  problem  description.  These  restrictions  are  made  only  to  simplify  the 
notation  and  the  formulation  of  the  equations,  but  the  scheme  can  be  easily  extended  to 
more  bodies.  Details  regarding  a  straightforward  implementation  to  a  system  of  multiple 
bodies  is  be  developed  in  Chapter  7. 

Figure  2.1,  is  a  schematic  drawing  of  two  elastic  bodies  in  their  reference  and  deformed 
configurations.  Each  body  is  represented  by  an  open  set  for  ?  =  1,2,  and  both  reside 
in  K  ,  where  is  the  number  of  space  dimensions.  These  bodies  undergo  motions 
described  at  time  t  by  the  deformations  for  i  =  1,2,  where  t  belongs  to  the  time 

interval  I  =  [0,  T].  Thus,  these  motions  are  expressed  via  the  following  mappings: 


:  12^  X  I  ^  E”"*™  for  z  =  1, 2  .  (2.1) 

For  simplicity,  we  assume  that  the  reference  configurations  12^*)  for  z  =  1,2  are  the  initial 
configurations  of  the  bodies  at  t  =  0.  We  also  assume  that  at  time  t  =  0,  there  is  no  contact 
between  the  bodies;  consequently,  no  contact  forces  are  present. 

We  identify  the  material  particles  of  each  solid  with  the  reference  coordinate  X  G  C 

.  The  current  placement  of  the  material  particle  X  €  12^^  at  time  t  G  [0,  T]  is  then 
expressed  as  :=  (p^^^{x,t). 

We  denote  the  boundaries  of  each  solid  i,  for  z  =  1, 2,  by  :=  512(*)  in  the  reference 
configuration  and  by  :=  in  the  current  configuration.  We  denote  by  %  := 

the  common  current  boundary  in  contact  between  any  two  bodies  (see  Figure  2.1); 
analogously,  in  the  reference  configuration  the  common  boundary  at  time  t  can  be  expressed 

byr,:=vr'(7a). 

We  denote  by  P<')  the  nominal  stress  tensor  (first  Piola-Kirchhoff  stress  tensor)  in  each  solid 
and  we  shall  restrict  our  analysis  to  hyperelastic  solids  characterized  by  their  respective 
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Figure  2.1;  Schematic  drawing  of  the  motion  of  two  bodies 

stored  energy  functions  where  =  Gradv?^*^  is  the  deformation  gradient  of  each 

solid.  Thus,  we  have 


pii)  = 


dW^i) 


(2.2) 


To  satisfy  the  principle  of  material  frame  indifference,  the  stored  energy  function  is 
invariant  under  the  action  of  the  proper  orthogonal  group  (the  rotation  group)  SO{ndim), 
that  is. 


1F«(QF«)  =  TT«(F«)  VQ  €  SOiudim) . 
Considering  a  one-parameter  group  of  rotations  Q{r])  with 


(2.3) 


—  W  G  Sofjldim)^ 


(2.4) 


77=0 
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where  so{n(iim)  denotes  a  linear  space  of  skew-symmetric  tensors.  Taking  the  derivative  of 
(2.3)  with  respect  to  77  and  setting  77  =  0,  we  obtain  the  following  relation: 

:W  =  0  yw  e  soiridim)  ,  (2.5) 

implying  the  symmetry  of  the  Kirchhoff  stress  tensor,  that  is: 


.^(0 ._  -T-m 

dF(^  ~ 


(2.6) 


The  symmetry  relation  2.6  leads  to  the  classical  conservation  law  of  the  total  angular  mo¬ 
mentum  as  discussed  below.  Furthermore,  a  classical  argument  (see  e.g.  Truesdell  & 
Noll[34])  leads  to  the  dependence  of  the  stored  energy  function  on  the  Green-Lagrange 
strain  tensor  defined  as 


£  = - 1) . 


(2.7) 


The  new  expression  for  the  stored  energy  function  denoted  by  W  is  then 


PF(*)(f(*>)  =  W(*)(f(*))  .  (2.8) 

As  an  example,  most  of  the  simulations  presented  in  this  work  use  the  Saint- Venant  Kirchhoflf 
model,  characterized  by: 


W{F)  =  W{E)  =  \e:CE,  (2.9) 

where  C  denotes  the  material  secant  tangent. 

We  denote  by  :=  the  material  velocity  field  of  the  solid  i,  and  the  corresponding 
reference  density.  The  superimposed  dot  (.)  refers  to  the  (material)  derivative  with  respect 
to  time  t. 
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2.3  Local  governing  equations 

In  the  absence  of  contact,  the  strong  form  of  the  local  momentum  balance,  along  with  the 
boundary  and  initial  conditions,  is  as  follows: 

For  all  t  6  [0,T],  [%  =  i,  2)  must  satisfy 


(*) 

pV 


pW 


^b)  =  Divp(*)  +  £)t  in 

(2.10) 

on  r^\ 

(2.11) 

on  r», 

(2.12) 

and  the  initial  conditions 


‘P 


(i) 


t=0 


dt 


<p 


(i) 


t=0 


1  (the  identity  mapping)  in 
in 


(2.13) 

(2.14) 


In  the  above  equations,  pQ^  is  the  reference  density,  tiq ^  is  the  outward  normal  in  the  reference 
configuration  and  represents  the  first  Piola-Kirchhoff  stress  tensor  for  body  z;  and  Div 
represents  the  divergence  operator  in  the  reference  geometry.  Also,  the  mappings  x 

I  — >  X  I  — >  MP'dim ^  ^(t)  .  pb)  x  I  represent  the  prescribed  body 

forces,  tractions  and  displacements,  respectively,  expressed  in  the  reference  configuration. 
Equation  2.13  is  the  mathematical  expression  which  states  that  the  initial  configuration  of 
the  bodies  coincides  with  the  reference  configuration.  In  equation  2.14,  represents  the 
prescribed  initial  velocity  field. 

For  the  above  system  of  equations  to  describe  a  well-posed  problem,  we  need  to  satisfy  the 
following  relations: 
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r^)ur^^ur^'>  =  (2.15) 

/^(i)nr»  =  r«n4*)  =  =  0.  (2.16) 

A  purely  Neumann  problem  has  no  displacements  imposed  on  the  boundary  of  its  bodies; 
thus  rjp  =0.  A  purely  Dirichlet  problem  has  no  tractions  imposed  on  the  boundary  of 
its  bodies;  thus  =  0.  In  this  dissertation,  we  work  exclusively  with  Neumann  type  of 
problems  as  they  have  very  distinct  conservation  properties  which  we  want  our  algorithms  to 
emulate.  However,  the  schemes  developed  in  this  dissertation  are  not  in  any  way  restricted 
to  this  type  of  problems. 

2.3.1  Contact  definition 

We  now  consider  the  addition  of  the  contact  constraints  to  the  strong  form  presented  in  the 
previous  section.  To  be  able  to  do  this  we  first  need  to  define  a  measure  of  the  closeness 
of  two  bodies.  This  is  done  by  means  of  the  gap  function.  The  gap  function  measures  the 
distance  between  one  point  on  the  surface  of  body  1  to  a  point  on  the  surface  of  body  2. 
This  mapping  between  a  point  of  body  1  and  and  its  closest  point  in  body  2  is  called  the 
closest  point  projection  and  its  definition  is  detailed  in  the  following  section.  The  concept  of 
closest  point  projection  with  its  algorithmic  interpretation  has  been  extensively  documented. 
Several  examples  can  be  seen  in  Hallquist  et  al[12]  and  Wriggers[36].  The  definition  of 
the  closest  point  projection  in  the  continuum  setting  is  documented  in  Laursen  &  Simo[19] 
and  in  SiMO  &:  Laursen[29]. 

2. 3. 1.1  Closest  point  projection 

In  what  follows,  we  denote  by  JC  a  material  point  which  belongs  to  the  surface  of  body  1, 
i.e.  to  r^.  We  define  the  gap  function  g{x)  for  a  material  point  X  e  F^  as  follows: 
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g{x)  :=  u .  -  ^^(y(x))]  ,  (2.17) 

where  the  mapping  Y  =  Y{x)  €  defines  the  closest  point  projection  on  the  contact 
surface  at  the  current  configuration  of  the  solids,  that  is, 

y(x)  =  argmin{||vi(x)-^2(y)|||  (2.18) 

yer2nrc 

In  the  above  equation,  ||  •  ||  denotes  the  usual  Euclidean  vector  norm  and  v  =  i/(y(x)) 
denotes  the  unit  outward  normal  to  the  current  contact  boundary  7^  D  7c. 

2.3.2  Contact  constraint 

The  normal  contact  constraint  enforces  the  physical  condition  of  impenetrability  and  com¬ 
pressive  interaction  between  bodies.  For  this  purpose,  we  define  one  of  the  surfaces  as  the 
slave  surface  or  contactor  surface  and  the  other  as  the  master  surface  or  target  surface.  We 
can  interpret  from  the  above  that  the  master  surface  nodes  (master  nodes)  define  the  surface 
of  impenetrability  for  any  node  belonging  to  the  slave  surface  (slave  node). 

For  simplicity,  we  always  choose  the  slave  surface  to  be  and  the  master  surface  to  be 
The  impenetrability  constraint  is  then  established  as 


g{X,t)>0  yxerU  (2.19) 

The  unilateral  constraint  is  intimately  related  to  the  definition  of  the  gap  function  g{X,t) 
through  the  closest  point  projection.  This  definition  is  fairly  straightforward  for  the  contin¬ 
uum  but  gives  rise  to  many  options  when  we  work  with  a  space  discretized  formulation  such 
as  the  finite  element  formulation.  One  of  these  options  is  developed  in  this  work.  Figure  2.2 
depicts  a  violation  of  condition  2.19  and  gives  a  physical  interpretation  of  the  definition  of 
the  gap  function  g{x,t). 
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Figure  2.2:  Geometric  interpretation  of  the  gap  function  g 

Remark  2.1  In  the  continuum  setting  the  distinction  between  master  and  slave  surfaces 
breaks  down,  because  if  one  point  belonging  to  the  slave  surface  lies  on  the  other  side  of 
the  master  surface,  the  reverse  is  automatically  true.  But  in  a  discretized  space  setting,  one 
can  find  many  cases  where  the  penetration  state  would  be  altered  if  we  interchanged  the 
definitions  of  slave  and  master  surface.  We  develop  this  idea  in  Chapter  4  when  we  explain 
the  finite  element  implementation  of  the  contact  problem. 

When  two  bodies  come  into  contact,  contact  forces  are  generated  along  the  common  bound¬ 
ary  7c.  We  recognize  two  contributions  to  this  traction:  a  normal  component  which 
originates  from  the  impenetrability  constraint  and  a  tangent  component  tj-  originating  from 
the  frictional  phenomena.  Thus,  we  have 


f(»)  :=  (2.20) 

=  4^1/  -  .  (2.21) 

To  enforce  the  principle  of  action  and  reation  across  the  the  contact  interface,  so  as  to 
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conserve  linear  momentum,  we  require  that  the  differential  contact  force  induced  on  the 
surface  of  body  2  at  F  is  equal  and  opposite  to  that  produced  on  body  1  at  X.  That  is, 


fdr^  +  i^dr^  =  0  ,  (2.22) 

where  is  the  traction  vector  due  to  contact  at  the  common  boundary  7c  belonging  to  the 
surface  of  each  body  i.  In  addition,  we  assume  that  no  tensile  (normal)  tractions  occur  due 
to  contact,  so  that 


>  0  ,  (2.23) 

where  is  the  unit  outward  normal  vector  belonging  to  a  point  on  the  surface  of  body  i. 
This  condition  ensures  that  the  two  surfaces  will  not  become  “glued”  together  once  contact 
occurs.  The  inequality  2.23  can  be  rewritten  as 


tjv  >0. 


(2.24) 


2.3.3  Persistent  contact 

Persistent  contact  occurs  when  two  surfaces  remain  in  contact  for  a  period  of  time,  i.e. 
contact  and  release  do  not  occur  instantaneously.  During  the  period  of  contact,  then,  all 
the  higher  time  derivatives  of  the  gap  function  g  vanish  (see  e.g.  Taylor  &  Papadopou- 
los[33]).  That  is, 


for  n  >  1  . 


In  particular,  one  can  define  a  gap  velocity  function  g  as 


(2.25) 
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g  =  [(p\x)  -  <p^{Y{x))]  ■  u  +  [¥>^(X)  -  <^^(r(x))]  •  i> .  (2.26) 

Since  differentiating  the  unit  normal  vector  u  in  time  leads  to  a  vector  which  is  perpendicular 
to  it,  we  then  have 

g{XN)  =  [i>\x)-,i>'^{Y{X))]-v 

=  [V^(X)  -  V^{Y{X))]  •  1/ .  (2.27) 

During  persistent  contact  the  gap  remains  constant  in  time,  i.e.  p  =  0;  thus  the  gap  velocity 
function  g  is  zero  during  this  period  of  time. 

Drawing  an  analogy  with  the  elasto-plastic  problem,  as  seen  in  Laursen  &  Simo[20],  we 
can  express  the  contact  problem  including  the  persistency  constraint  in  the  following  way: 


9{X,t) 

>  0, 

(2.28) 

tN{X,t) 

>  0, 

(2.29) 

tNix,t)gix,t) 

=  0, 

(2.30) 

tN{X,t)9{X,t) 

=  0. 

(2.31) 

Equations  2.28-2.30  are  the  Kuhn-Tucker  complementary  conditions  and  equation  2.31  is  the 
persistency  condition.  Thus,  equations  2.28-2.30  reflect  the  impenetrability  constraint,  the 
compressive  normal  traction  constraint  and  the  requirement  that  the  pressure  be  nonzero 
only  when  contact  occurs  (i.e.  when  ^  =  0).  The  persistency  condition  states  then  that  the 
rate  of  separation  between  the  two  surfaces  will  be  nonzero  only  when  the  contact  pressure 
vanishes. 
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2.4  The  frictional  problem 

Throughout  this  work,  we  consider  the  Coulomb  friction  model  with  no  evolution  or  rate 
dependence  of  the  coefficient  of  friction  fi,  although  such  effects  may  be  included  in  the 
presented  methodology  with  slight  modifications  (see  e.g.  Wriggers  &  Stein[37]). 

2.4.1  Contact  kinematics 

A  mathematical  description  of  the  relative  motion  between  two  bodies  in  contact  is  essen¬ 
tial  for  the  consistent  statement  of  the  law  governing  the  friction  phenomena.  Below,  we 
shall  describe  the  necessary  components  of  contact  kinematics,  also  given  in  Wriggers  & 
Stein[37]  and  Laursen  &:  Simo[19].  Through  the  definition  of  the  gap  function  g  by  means 
of  the  closest  point  projection,  we  define  a  convected  basis  along  the  master  surface  which 
enables  us  to  express  the  frictional  constraints  Parametrizations  for  and  7^  are  shown  in 
Figure  2.3. 


Figure  2.3:  Parametrization  of  the  contact  surfaces  and  7^. 

The  convected  bases  for  and  7^,  that  is,  the  tangent  vectors  in  the  reference  and  current 
configuration,  are  defined  via  partial  derivatives  with  respect  to  the  parametrization  variable 
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^  in  the  parent  domain  as 


and 


Ta{^)  ■■=  Ka(^) 


(2.32) 


r-a(4)  ■■=  =  F^(^o"(^))j^aU)  ,  (2.33) 

respectively,  for  cc  =  1,2,  where  is  the  deformation  gradient  corresponding  to  the  defor¬ 
mation  (p^(x).  In  equations  2.32  and  2.33  the  expression  (•)^q  means  the  partial  derivative 
with  respect  to  ^q,  where  a  =  1,2.  In  a  two  dimensional  setting,  the  tangent  plane  is  one 
dimensional,  hence  a  is  omitted. 

For  any  point  X  G  we  obtain  a  corresponding  point  Y  e  through  the  closest  point 
projection  minimization  indicated  in  equation  2.18,  that  is  y  =  Y'(x).  This  calculated 
contact  point  is  expressed  in  a  parametrized  form  for  the  reference  and  current  configurations 
as 


Y{x,t)  =  ^^{i{X,t))  (2.34) 

and 

y{X,t)=^Ui{X,t))  ,  (2.35) 


respectively,  for  o;  =  1, 2. 

We  denote  by  Ma^  and  the  components  of  the  associated  positive  definite  metrics  for 
the  reference  and  current  configurations,  respectively,  and  their  expressions  are  given  by 


(2.36) 
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and 


:=  Ta-rp  ,  (2.37) 

for  a,  ^  =  1, 2.  As  mentioned  above,  for  a  two  dimensional  contact  problem  these  metrics 
become  positive  scalars. 

The  friction  phenomenon  is  ruled  by  the  amount  of  relative  motion  between  the  surfaces. 
Thus,  to  measure  the  amount  of  relative  slip,  we  need  to  calculate  the  time  derivative  of  the 
material  contact  point  y  €  the  relative  slip  velocity,  obtained  through  the  closest  point 
projection,  i.e.. 


Y  =  ,  (2.38) 

expressed  in  terms  of  the  basis  in  the  reference  configuration.  By  taking  the  time  derivative 
of  the  expression 


ip\x)  -  v^{Y{x))  =  gv  ,  (2.39) 

we  obtain 

Ivj  -F'^Y  =  gu-g  ^  ^2.40) 

where  the  relative  material  velocity  is  expressed  as 

|yl  :=  yi(x)  -  v'^{Y{X))  ,  (2.41) 

with  :=  for  z  =  1, 2.  In  general,  we  use  !•]]  to  express  the  relative  difference  of 
the  variable  (•)  evaluated  at  the  two  contact  points  X  and  Y.  The  normal  and  tangent 
components  of  expression  2.40  are,  respectively, 


9  =lVl-v 


(2.42) 
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and 

=  1^1  -Ta  +  gV^-i^,  (2.43) 

where  the  symmetric  matrix  A  has  the  following  expression: 

Aai3  :=  mai3  -  gTa,i3  ■  v  .  (2.44) 

The  matrix  A  is  assumed  to  be  invertible  at  all  times  so  that  equation  2.43  defines  the  slip 
rate  uniquely.  Notice  that,  in  the  limit  when  ^  =  0,  this  property  follows  immediately 
from  the  positive  definiteness  of  the  metric  ma^.  Since  it  is  this  particular  limit  case  which  we 
are  trying  to  enforce  throughout  the  problem,  the  previous  assumption  is  not  too  restrictive. 

2.4.2  Coulomb  friction 

The  friction  contribution  to  the  surface  traction  due  to  contact  is  expressed  as  follows: 


tx  =  —txpT^  ■  (2.45) 

Coulomb  friction  is  described  by  means  of  a  slip  surface  defined  by  the  following  slip  function 


(f>:=\\tT\\-pitN .  (2.46) 

The  behaviour,  under  Coulomb  friction,  is  characterized  by  a  perfect  stick  condition  (static 
friction)  until  the  value  of  the  slip  function  <j)  becomes  zero  due  to  the  increase  of  the  tan¬ 
gential  traction  or  to  the  decrease  of  the  normal  traction  at  that  point.  Thereafter,  frictional 
slip  (dynamic  friction)  occurs  where  the  slip  is  related  to  the  tangential  traction  (see  e.g. 
Laursen  &  Siiv[o[20]),  as  follows: 


VT  = 


(2.47) 
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The  slip  can  also  be  defined  by  means  of  kinematic  considerations  (see  e.g.  Laursen  & 
Simo[20]),  as  follows: 

:=  f^iAnx))) 

=  .  (2.48) 

The  frictional  behaviour  can  then  be  described  using  Kuhn-Tucker  complementary  conditions 
and  the  persistency  condition: 


<l>  =  ||tr||  -  jutiv  <  0  , 

(2.49) 

IV 

o 

(2.50) 

o 

II 

(2.51) 

II 

o 

(2.52) 

where  the  Euclidean  norm  in  is  given  by  HirlP  =  m°‘HTatTp  Laursen  &  Simo[20]. 

Remark  2.2  As  recognized  in  Michalowsky[22],  the  slip  law  expressed  in  equation  2.47 
is  non-associative  due  because  the  slip  potential,  whose  gradient  defines  the  slip  direction,  is 
not  the  same  as  the  slip  function  (f).  If  these  were  the  same,  the  slip  would  have  a  dilatational 
component  due  to  the  dependency  of  (f)  on  the  normal  contact  pressure  tjv. 

In  order  to  approximate  the  perfect  stick  condition  until  the  tangential  traction  reaches 
a  certain  value  proportional  to  the  normal  pressure,  we  perform  a  penalty  regularization 
of  equation  2.47  (see  e.g.  Wriggers  &  Stein  [37]  and  Laursen  &  Simo[20]).  This 
regularization  is  achieved  through  the  introduction  of  the  time  rate  of  change  of  the  frictional 
traction,  in  a  basis  which  must  be  suitably  defined. 

Frame  indifiference  is  achieved  through  the  choice  of  the  convected  basis  in  which  all  expres¬ 
sions  are  written.  This  property  will  continue  to  hold  if  the  rate  of  change  is  expressed  using 
the  spatial  Lie  derivative  of  the  frictional  traction: 
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CtT  :=  .  (2.53) 

We  introduce  the  (large)  tangential  penalty  parameter  kt  for  the  regularization  of  the  fric¬ 
tional  response  as  follows: 


^  <  0  , 

(2.64) 

tT  1  ^ 

=  7 11,  11+  CtT, 

1I«t||  Kt 

(2.55) 

IV 

o 

(2.56) 

II 

o 

(2.57) 

Using  the  Lie  derivative  in  the  regularization  permits  the  system  of  equations  2.54-2.57  to  be 
solved  in  closed  form  as  in  the  elastoplastic  problem,  with  the  tangential  penalty  parameter 
kt  playing  the  role  of  the  elastic  modulus.  The  above  statement  of  the  problem  allows  some 
slip  to  occur  before  frictional  slip  initiates.  The  amount  of  premature  slip  depends  on  the 
magnitude  of  the  tangent  penalty  parameter  kt  and  will  vanish  as  kt  oo.  One  can  make 
two  points  regadrding  the  analogy  between  frictional  and  elastoplastic  response: 

1.  The  tangential  penalty  is  used  to  preclude  elastic  tangential  displacement,  while  in 
elastoplasticity  it  is  used  to  regularize  rigid-plastic  response. 

2.  It  has  been  suggested  that  the  tangential  stiffness  kj-  is  related  somehow  to  the  elastic 
stiffness  of  contacting  asperities,  so  that  there  is  a  physical  meaning  to  the  analogy 
with  an  elastic  modulus. 
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Chapter  3 


The  Weak  Formulation 


3.1  Introduction 

In  this  chapter,  the  local  equations  presented  in  Chapter  2  are  incorporated  into  the  global 
variational  formulation.  As  in  Laursen  &  Simo[20],  the  weak  form  of  the  initial  boundary 
value  problem  includes  the  local  governing  equations  and  their  associated  kinematic  quanti¬ 
ties  directly  so  that  one  obtains  a  variational  principle  for  the  two  body  problem  that  is  a 
geometrically  exact  statement  of  the  contact  virtual  work. 

We  also  make  use  of  the  global  statement  to  investigate  the  properties  of  the  underlying 
continuum  problem,  by  studying  the  behaviour  of  quantities  such  as  linear  momenta,  an¬ 
gular  momentum  and  total  energy  of  the  system.  This  provides  us  with  the  guidelines  to 
later  develop  the  algorithmic  counterpart  of  the  weak  form  so  that  it  inherits  its  physical 
properties. 

3.2  Weak  form  of  the  governing  equations 

The  weak  form  of  the  momentum  balance  equations  can  be  expressed  as  follows: 
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V'  f  pV)y{^  dO  +  y  f  :  Grad(<5v^(*))  dn 

i  Jnii)  i  Jnii) 

=?L  •  5<p^^  dn  +  i-Scp^^  dr 

+  f  [tN^g  — 

Jure 


for  all  admissible  variations  5ipV)  (i=  2)  satisfying 


The  terms  on  the  right-hand  side  of  equation  3.1  represent  the  work  done  by  the  body  forces 
b,  the  surface  tractions  i  and  the  contact  tractions  i. 

By  manipulating  the  definition  of  the  gap  function  g  as  expressed  in  equation  2.17,  we 
obtain  the  expression  for  5g,  so  that  by  replacing  the  time  derivative  in  equation  2.42  with 
the  variation  of  we  obtain 


7] :  ->■  I”''*” 


^1 


r>»  =  ®} 


(3.2) 


Sg  =  [<^1  •  i' 

=  ^  .  (3.3) 

The  same  replacement  can  be  done  in  equation  2.43  to  find  5^“  for  a  =  1, 2,  yielding  the 
following  expression: 


. 

=  -  (5¥>^(y(X))]  .  -h  g5^\  •  I.  ,  (3.4) 

where  Aag  has  the  same  definition  as  in  equation  2.44.  Notice  that  during  contact,  the 
solution  has  to  satisfy  p  =  0,  hence  we  can  rewrite  equation  3.4  as 
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-  S<p‘^{Y{x))]  •  .  (3.5) 

So  far  we  have  followed  closely  the  geometric  structure  of  the  contact  constraints  and  have 
not  made  any  assumptions  which  are  not  valid  for  the  continuum. 


3.3  Properties  of  the  weak  formulation 

In  this  section,  we  study  the  evolution  of  various  physical  quantities  associated  with  the 
mechanical  system.  In  particular,  we  observe  that  the  frictionless  contact  problem  among 
elastic  bodies  allowing  finite  deformations  is  a  unilaterally  constrained  infinite  dimensional 
Hamiltonian  system  which  leads  to  the  conservation  laws  described  below. 

Consider  the  total  linear  momentum, 


L  :='y  f  , 


and  the  total  angular  momentum. 


J  := 


vjW  X  p<f'v<'‘'>dn , 


(3.6) 


(3.7) 


of  a  pair  of  solids  (i  =  1,2).  The  symbol  x  refers  to  the  cross  product  of  two  vectors  in 
or  the  equivalent  definition  if  we  work  in  a  two  dimensional  formulation.  Similarly,  we 
consider  the  total  energy  of  the  system. 


e 


.  v^^dQ  +  y^  f  dQ 


)C  +  W, 


(3.8) 


where  the  total  kinetic  energy  of  the  system  is  denoted  by  K  and  the  total  elastic  strain 
energy  of  the  system  is  denoted  by  W. 
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The  case  of  interest  in  this  analysis  is  the  homogeneous  Neumann  problem,  for  which  neither 
boundary  displacements  nor  external  loading  is  imposed;  consequently,  the  linear  and  angular 
momenta  are  conserved  (i.e.  L  =  constant  and  J  =  constant). 

For  the  frictionless  problem,  tx  =  0,  the  linear  momentum  L,  the  angular  momentum  J  and 
the  total  energy  of  the  system  S  are  conserved  (i.e.  these  quantities  are  constants  of  motion). 
When  friction  is  present,  energy  dissipates  from  the  system.  These  conservation  properties 
are  investigated  in  the  next  section. 

3.3.1  Conservation  properties  of  contact 

3. 3. 1.1  Conservation  of  linear  momentum 

Since  no  boundary  displacements  are  imposed  in  this  case,  =  0  and  we  can  use  a  constant 
as  an  admissible  variation,  as  follows: 


for  «  =  1,2,  with  a  €  constant.  For  this  case,  then,  Grad(5¥5^*^)  =  0.  Using  equation 
3.1  together  with  the  definition  3.6  and  taking  into  account  that  f  =  0  and  6  =  0,  we  can 
state  the  following: 


•  a  dQ 
a]  dP 


.0  ,,yfl  €  . 


(3.10) 


Consequently,  ^  =  0. 


_ _ 3.3.  PROPERTIES  OF  THE  WEAK  FORMULATION 

3.3.1. 2  Conservation  of  angular  momentum 
Let  us  consider  the  following  admissible  variations: 

(3.11) 

for  i  =  l,2,  with  w  €  E"*™  constant,  and  We  then  obtain 

Grad(5yj®)  =  ,  (3.12) 

where  W  is  the  skew  symmetric  tensor  cooresponding  to  the  axial  vector  to,  that  is, 


Wa  =  wXa  Va  €  .  (3.13) 

Substituting  the  variation  defined  in  equation  3.11  into  equation  3.1,  and  using  the  definition 
3.7  and  the  equality  3.12,  we  obtain  the  following  expression: 


=  Y^w-  [  X  X 

=  dQ 

=  Y  f  •  (to  X  x(*))  dn . 


(3.14) 


Also,  for  this  particular  variation  Sg  takes  the  following  form: 


5g  =  1/  •  [to  X  {,p^{x)  -  (y(A:)))] 

=  gu  •  [w  X  i/] 

=  0  Vto  G  .  (3.15) 
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where  we  have  used  the  definition  of  the  gap  function  p  in  2.17.  For  this  choice  of  variation, 
the  expression  6^^  takes  the  following  form: 

=  To,  -  [w  X  {<p^{x)  -  tp^  (y(x)))]  Pgv  [to  X 

=  9To,  -[w  X  1>]  +  gu-[w  X  Ta] 

=  0  Vto  e  R"**'”*  ,  (3.16) 

where  we  have  once  more  used  the  definition  of  the  gap  function  g  in  equation  2.17.  Using 
these  last  two  results,  i.e.  equations  3.15  and  3.16,  we  can  express  the  variation  of  the  angular 
momentum  as  follows: 


-  /  dQ 

J  u/2(0 

0  Vto  €  ,  (3.17) 

where  we  have  used  the  symmetry  property  of  the  Kirchhoff  stress  tensor,  detailed  in  equation 
2.6.  Consequently,  ^  =  0. 

3.3. 1.3  Evolution  of  the  energy 

The  evolution  of  the  energy  can  be  obtained  by  choosing  the  following  variation: 


for  i  =  1, 2  . 


(3.18) 


Substituting  this  variation  into  equation  3.1  and  using  the  definition  3.8,  we  obtain  the 
following  expression: 
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dt 


Jure 

-  f  trgi^  dr  , 

Jure 


liW  '■  Gtad(l^«)  dn 


(3.19) 


where  we  have  used  the  definition  of  the  gap  velocity  function  g  expressed  in  equation  2.42 
and  the  persistency  condition  referred  to  in  equation  2.31.  From  expression  3.19,  we  can 
observe  that  any  change  in  the  energy  of  the  system  comes  from  the  frictional  term.  Energy 
is  dissipated  only  when  frictional  slip  occurs,  because  of  the  condition  7(?!)  =  0  in  2.57.  That 
is,  if  (^  <  0,  then,  7  =  0.  Using  this  last  result  in  equations  2.47  and  2.48,  we  find  that  ^  =  0. 
When  frictional  slip  occurs,  0  =  0  due  to  condition  2.57;  therefore,  combining  equations  2.47 
and  2.48  yields  an  expression  for  £,  given  by 


INI 


(3.20) 


Substituting  equation  3.20  into  equation  3.19,  and  observing  that  Utrll  =  we  find 


<  0, 


where 


(3.21) 


ll*r||re/  •=  tTaM°‘^tTg  ■ 


(3.22) 


We  have  used  equations  2.24  and  2.56  to  establish  the  sign  of  the  integrand. 

We  denote  by  Vfric,  the  energy  dissipation  introduced  by  the  friction  phenomena: 


>  0. 


(3.23) 
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Consequently,  §  =  -Vfric- 

Remark  3.1  Notice  that  in  the  absence  of  friction,  namely  //  =  0,  the  energy  of  the  system 
is  conserved  (i.e.  ^  =  0). 
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Chapter  4 


Finite  Element  Implementation 


4.1  Introduction 

The  finite  element  method  is  applied  to  discretize  the  equations  delineated  in  the  previous 
chapters,  including  the  contribution  of  the  contact  forces,  though  the  actual  expression  of 
the  normal  contact  pressure  and  tangential  traction  is  posponed  until  Chapters  5  and  6). 

The  finite  element  implementation  developed  in  this  chapter  is  equivalent  to  those  presented 
by  Wriggers  &  SiMO  [36]  for  the  two  dimensional  case  and  by  Parisch[27]  for  the  three 
dimensional  case.  This  section  is  based  particularly  on  the  work  done  by  Laursen  and 
SiMO  [20],  who  developed  a  general  methodology  for  isoparametric  elements  in  two  and 
three  dimensions. 


4.2  Time  discretization  of  the  weak  equation 

We  perform  a  discretization  of  the  time  interval  of  interest  [0,  oo]  into  subintervals  [t„,  tn+i]- 
We  denote  by  At  >  0,  the  corresponding  time  step  At  =  tn+i  —  tn,  where  we  assume  that 
tn  =  nAt.  For  a  variable  g(t)  continuous  in  time,  we  perform  an  algorithmic  approximation: 
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Qn  «  g{nAt)  =  g{tn) 


(4.1) 


and 


3n+^  2  9n+l)  • 


(4.2) 


In  the  absence  of  external  forces  and  imposed  displacements,  equation  3.1  is  discretized  in 
time  as  follows: 


y  f  pn+Ui)  ,  Grad((5^(*))  dQ 
=  f  [tj^Sg  -  dr  ,  (4.3) 

JVFr. 


for  all  admissible  variations  i  =  l,2. 


Following  the  work  of  SiMO  &  Tarnow  [30],  the  internal  elastic  forces  can  be  time  dis¬ 
cretized  in  such  a  way  as  to  inherit  all  the  conserving  properties  inherent  to  the  unconstrained 
continuum  problem.  The  algorithm,  which  allows  finite  deformations,  is  a  second  order 
conserving  approximation  of  the  internal  force  term  in  equation  3.1  for  the  elastodynamic 
problem,  and  is  stated  as  follows: 


P'1+1. (i) 


Grad(5v’^*^)  dH  =  f  V„+i  (^^p)  :  dQ  , 


(4.4) 


where  the  Kirchhoff  stress  r  at  time  t„_,.  i  is  expressed  using  using  the  Saint- Venant  Kirchhoff 
model  detailed  in  equation  2.9: 


^(n+i)  _  {lc{En  +  £;„+i))  .  (4.5) 

The  deformation  gradient  tensor  is  computed  at  the  mid-point  configuration  1 ,  and 

2  2 

is  given  as  follows: 
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:=  Grad^„+i  .  (4.6) 

The  Green-Lagrange  strain  tensors  En+i  and  En  are  computed  at  the  configurations  corre¬ 
sponding  to  tn+i  and  tn,  respectively. 

Thus,  we  can  rewrite  equation  4.3  as 


f  :  lc{En+i  +  En)  dQ 


=  I  ti<^5g  —  dr  , 
Jure 


(4.7) 


combined  with  the  following  standard  time  stepping  scheme: 

-  <Pn)  =  . 


(4.8) 


4.3  Conserving  properties  of  the  algorithm 


In  SiMO  &  Tarnow  [30],  this  algorithm  is  shown  to  have  important  conserving  properties, 
which  we  detail  below.  In  the  following  proofs,  we  assume  that  we  have  a  homogeneous 
Neumann  problem,  i.e.  that  no  displacements  are  imposed  on  the  body  (jT^  =  0)  and  no 
external  body  forces  or  surface  tractions  are  applied  to  the  body  (6=0  and  t  =  0).  Most 
importantly,  we  assume  that  the  problem  is  unconstrained,  that  is  that  no  contact  forces  are 
present  as  yet  in  the  problem  {Pc  —  0). 


4.3.1  Internal  linear  momentum 

In  the  time  discretized  setting,  conservation  of  linear  momentum  holds  under  the  above 
assumptions  and  is  expressed  in  the  following  way: 
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^n+l  —  ^ 


where  the  total  linear  momentum  is  defined  by 


Lt  := 


(4.9) 


(4.10) 


Using  the  expression  4.7  and  choosing  =  a  e  where  a  is  any  constant, 


[•^n+l  Lji  ■  O, 


(4.11) 


since  (o)  =  0.  We  conclude  that  there  is  no  contribution  of  the  internal  forces  to  the 
linear  momentum  Lt. 


4.3.2  Internal  angular  momentum 

In  the  time  discretized  setting,  conservation  of  angular  momentum  holds  under  the  above 
assumptions  and  is  expressed  in  the  following  way: 


•7n+l  —  Jfi  ) 

where  the  total  angular  momentum  of  the  system  is  defined  by 


(4.12) 


We  can  reorganize  the  previous  expression  in  the  following  way: 


(4.13) 
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Ln+l  Tn  — 


V’Si  X  -  ¥>«  X 


dQ 


=  -V®) 


+ 


n+i 


dQ  . 


(4.14) 


Using  equation  4.8,  we  see  that  on  the  last  term  of  the  integrand  of  equation  4.14  vanishes 
due  to  co-linearity.  We  also  use  the  following  admissible  variation: 


^  >  (4-15) 

with  w  G  (constant)  and  substitute  it  into  equation  4.7.  Then,  the  first  term  of  the 

integrand  of  equation  4.14  becomes 


[*^n+l  *7n.  ■  W 


=  E”  ■  X  p®  (v®,  -  V®)  dn 

=  J Grad(a:„+i  x  dQ 

=  -'yw  [  W:t("+2)  dQ 

,  Ja(i) 


=  0, 


(4.16) 


where  we  have  used  the  result 

Grad  x  w)  =  WF^^  ,  (4.17) 

in  which  W  is  the  skew-symmetric  tensor  with  axial  vector  w.  As  the  contraction  of  a 
symmetric  and  a  skew-symmetric  tensor  vanishes,  we  conclude  that  there  are  no  contributions 
of  the  internal  forces  to  the  total  angular  momentum  of  the  system. 


34 


4.3.  CONSERVING  PROPERTIES  OF  THE  ALGORITHM 


4.3.3  Internal  energy 

In  the  time  discretized  setting,  we  prove  that  energy  is  conserved  under  the  above  assump¬ 
tions.  The  energy  of  the  system  consists  of  the  kinetic  energy  /C*  for  t  €  [tn,in+i],  which  can 
be  written  as 


/Ct  = 


(4.18) 


and  the  potential  energy  Wt  for  t  e  caused  by  the  internal  elastic  forces,  given  by 


Wt  =  'y  f  W  dO, 


(4.19) 


where  W  is  the  elastic  potential  as  defined  in  equation  2.9.  Using  the  admissible  variation 


d(pV)  — 


(4.20) 


we  can  rewrite  the  left-hand  side  of  equation  4.7  in  the  following  way: 


=  [^n+i  -  ICn]  .  (4.21) 


where  ||  •  ||  denotes  the  standard  Euclidean  norm.  The  right-hand  side  of  equation  4.7,  when 
using  the  above  variation,  becomes 


=  -E/ 

V„+i  {cpn+\ 

2 

=  -e/ 

=  -(Wn+i- 

W„), 

1 

2^ 


(4.22) 
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where  we  have  used  the  Saint- Venant  Kirchhoff  model  to  describe  the  elastic  behaviour  of 
the  bodies.  We  conclude  then  that  the  total  energy  St  =  K,t  +  Wt  is  conserved. 

4.4  The  semi-discrete  equations 

The  weak  form  of  the  governing  equations  of  our  mechanical  system  is  discretized  in  space, 
resulting  in  a  set  of  semi-discrete  equations  as  explained  in  [14].  These  equations  are,  essen¬ 
tially,  a  set  of  nonlinear  ordinary  differential  equations. 

The  space  is  discretized  using  a  standard  isoparametric  finite  element  formulation;  thus 
equation  3.1  is  transformed  into  the  form 


d(t)  =  M-V(t)  , 

(4.23) 

P{t)  =  -fintid{t))  +  fc{d{t))  +  fext, 

(4.24) 

where  d  is  the  vector  of  nodal  displacements,  fint  is  the  vector  of  internal  forces  due  to  elastic 
deformation,  fext  represents  the  discretized  vector  of  external  forces  such  as  body  forces  and 
fc  is  the  vector  of  contact  forces.  We  have  introduced  the  vector  of  nodal  linear  momenta  as 
an  intermediate  quantity  which  is  defined  in  the  following  way: 


p  :=Mv  ,  (4.25) 

where  v  :=  d{t). 

As  for  any  initial  value  problem  the  equations  4.23  and  4.24  are  solved  with  a  set  of  initial 
conditions  for  d  and  d  (or  p),  given  by  .  .  ■  ^  ^ 


d(0)  =  do  , 
d(0)  =  Vo  ■ 
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(4.26) 

(4.27) 


4.5.  SPACE  DISCRETIZATION 


4.5  Space  discretization 

To  obtain  the  semi-discrete  system  of  equations  described  in  the  previous  section,  let  us 
consider  the  standard  isoparametric  finite  element  space  discretization, 


^node 

X=Y1  (4.28) 

A=1 

and 

Tlnode 

=  X+J2  >  (4-29) 

>1=1 

where  X  represents  the  reference  coordinate  field  and  represents  the  current  coor¬ 

dinate  field.  The  interpolation  shape  functions  are  expressed  by  N-^  :  □  — R  with  4  €  □ 
for  A  =  !,•••  ,nnode  (the  total  number  of  nodes  in  the  isoparametric  element).  Finally, 
we  have  the  vector  of  nodal  displacements  €  R”*™  and  the  vector  of  reference  nodal 
coordinates  Xa  €  .  The  variables  can  be  grouped  into  the  vector  d  €  R"'«  where 

Tleg  =  ^  '^node- 

The  mass  matrix  M  is  defined  by  the  standard  assembly  procedure; 


net 

M  =  (4.30) 

e=l 

where  M®  is  the  elemental  mass  matrix  and  is  the  total  number  of  elements.  For  an 
element  with  Uen  nodes,  the  elemental  mass  matrix  M®  is  as  follows; 


Mnindim 

^lUen  "^ndim 

(4.31) 


where  Inum  is  the  rank-two  identity  matrix  in  R"*'"  and  Mab  for  a  body  i  is  defined  as 
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Mab=  f  dO  (4.32) 

where  is  the  elemental  volume  for  body  i.  This  definition  of  the  mass  matrix  is  called 
the  consistent  mass  matrix.  In  subsequent  chapters  of  this  work,  we  may  also  use  the  lumped 
mass  matrix,  which  is  obtained  by  performing  a  standard  row-sum  technique  [14],  and  is 
expressed  as  follows: 


Mab  =  MaSab  (no  sum)  ,  (4.33) 

where  Ma  =  dQ. 

The  external  force  vector  f^xt  ^  1^"®'’  corresponds  to  the  contributions  of  the  body  forces  h 
and  the  imposed  external  tractions  t  applied  to  Tc,.  The  internal  force  vector  fint  € 
corresponds  to  the  stress  divergence  term  in  the  continuum;  its  expression  is  as  follows: 


fint  —  [ 

Junii) 

where  Bt  is  the  standard  strain  operator, 

BtSd  :=  =  sym  [GradfJte]  ,  (4.35) 

and  tt  :=  ¥>t(x)  —  JC  represents  the  displacement  field  at  time  t. 

4.5.1  Closest  point  projection  in  the  space-discretized  setting 

To  be  able  to  define  the  contact  force  vector  fc  in  equation  4.24,  we  first  need  to  develop 
the  concept  of  the  closest  point  projection  in  a  space  discretized  setting  and  expand  on  the 
notation  involved. 

We  make  use  of  the  standard  slave/master  data  structure  (see  [12]).  We  can  identify  a 
slave  body  and  a  master  body  whose  surfaces  are  now  composed  of  surface  elements  and 


38 


4.5.  SPACE  DISCRETIZATION 


surface  nodes.  We  perform  a  contact  detection  procedure  for  every  surface  node  in  the  slave 
body,  denoted  by  S,  to  check  for  possible  penetration  into  the  master  body.  A  contact  point 
Y  (Xs)  is  established  on  the  surface  of  the  master  body  which  can  be  identified  by  a  location 

within  a  particular  surface  element  defined  by  its  nodes,  i.e.  Ml,  M2, _ Figure  4.1  shows 

a  schematic  view  of  a  slave  node  S  penetrating  a  master  surface  element  and  the  definition 
of  the  contact  point  y(.X's)  by  means  of  a  closest  point  projection. 


Figure  4.1:  Schematic  drawing  of  the  closest  point  procedure  on  a  discrete  space  setting  in 
two  dimensions.  In  the  case  shown,  contact  is  detected  and  the  master  surface  element  is 
defined  by  4  nodes. 

The  closest  point  projection  shown  in  Figure  4.1  leads  us  to  establish  a  contact  element 
consisting  of  the  contacting  slave  node  S  (with  coordinates  Xs)  and  the  nodes  that  define 
the  master  segment  containing  the  penetration  point  Thus,  the  contact  element  can 

be  described  by  the  set  of  nodes  {S,  Ml,  M2,  M3,  M4, ....}. 

Remark  4.1  Throughout  the  evolution  of  the  problem,  the  number  of  contact  elements  will 
vary  depending  on  the  number  of  slave  nodes  contacting  the  master  body  at  any  given  time. 
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The  contact  force  vector  R  can  then  be  expressed  as  the  assembly  of  the  force  contributions 
from  each  contact  element,  as  follows: 


'^slave 

fc=  j\  fs,c  (4.36) 

S=rl 

where  risiave  is  the  number  of  contact  elements  (or  slave  node/master  segment  pairs)  at  any 
given  time  and  is  the  contact  force  applied  to  each  of  the  nodes  belonging  to  the  contact 
element.  It  is  defined  as 


Is,c  —  tj>fGg 


(4.37) 


where 


(4.38) 


and  i®  fhs  number  of  nodes  in  the  master  segment  in  contact  with  the  slave  node  S. 

In  equation  4.38,  N^^{^s)  denotes  the  standard  isoparametric  shape  function  of  node  MI  in 
the  master  segment  at  the  point  of  contact  with  normal  i/^,  obtained  by  the  closest  point 
projection  procedure;  see  Figure  4.1. 

The  discrete  counterpart  of  equation  2.39  is  now 


77  ^ 

**'master 

giXs)^  =  xs  —  )  (4.39) 

7=1 

where  xs  =  (p^{Xs)  and  xmi  =  <p‘^(Xs)  are  the  current  positions  of  the  slave  and  master 
segment  nodes,  respectively.  We  note  that  the  master  segment  shape  functions  satisfy  the 
relation 
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T>^ 

‘‘'master 

Y,  Ar"'(«)  =  l  (4.40) 

/=1 

for  any  point  ^  belonging  to  the  master  segment.  We  introduce  the  following  simplifying 
notation: 

ds  Vs 

dui  ,  ^  vm\ 

and  V5  = 

dM2  Vm2 

these  refer  to  the  displacement  and  velocity  of  the  slave/master  segment  nodes,  respectively. 
In  general,  we  denote  with  (.)j  the  quantities  that  refer  only  to  the  contact  element. 

Remark  4.2  As  mentioned  in  Chapter  2,  contact  detection  in  a  space  discretized  setting  by 
means  of  a  closest  point  projection  may  yield  very  different  results  Hallquist  et  al[12]. 
Figure  4.2  shows  an  example  where  choosing  which  surface  is  the  slave  reverses  the  contact 
state.  The  closest  point  projection  will  detect  contact  for  the  situation  depicted  on  the  left 
and  will  not  detect  contact  if  we  reverse  the  choice  of  slave  and  master  surface  as  shown  in 
the  figure  on  the  right. 

Me 

Figure  4.2:  Schematic  drawing  of  two  contact  situations  involving  with  two  dimensional 
linear  surface  elements.  Slave  nodes  are  denoted  by  S  and  master  nodes  by  M. 

4.6  Temporal  and  space  discretization 

We  perform  a  discretization  of  the  time  interval  I  into  a  series  of  time-steps,  i.e.  U^_o[tn)  in+i]- 
The  time  continuity  implied  in  expressions  4.23  and  4.24  is  removed  by  solving  the  ordinary 


Ma  Mb  Mb 
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differential  equations  via  a  time-stepping  algorithm  in  each  subinterval  [tn,  tn+i]-  We  consider 
a  mid-point  approximation  to  the  system  described  by  equations  4.23  and  4.24,  which  results 
in  the  following  system: 


^^(rfn+l  dfi)  —  I’n+i  > 

—M{vn+1  -  Vn)  =  -fint  "  +  /c  ^  +  f ext  , 

where  At  =  —  tn,  the  time  interval  of  interest.  We  have  used  the  following  approxi¬ 

mations:  dn  ~  d{tn),  Vn  ~  v(tn)  and  =  (vn+i  -l-v„)/2.  We  have  also  eliminated  the 
momentum  from  equation  4.23  by  using  the  following  approximation: 


(4.42) 

(4.43) 


p(t)  ^Pt  =  Mvt  for  t  e  Un[tn,  tn+l]  •  (4.44) 

The  discrete  form  of  the  normal  and  frictional  contact  forces  fc  are  discussed  and  expanded 
upon  in  the  next  chapter.  The  internal  force  vector  denoted  by  in  equation  4.43 

corresponds  to  the  time  discretization  presented  in  SiMO  &  Tarnow  [30].  This  algorithm 
is  a  second  order  conserving  approximation  of  the  internal  force  vector  fint  at  for  the 
elastodynamic  problem  which  allows  finite  deformations.  The  strain  operator  Bt  defined  in 
equation  4.35  is  evaluated  at  the  mid-point  configuration  v’n+l  ivn+i  +  v’n)/2  and  the 
internal  force  vector  is  then  given  by 


J int 


dfl  , 


(4.45) 


where  the  discrete  Kirchhoflf  stress  is  calculated  using  the  Saint- Venant  Kirchhoff 

model  defined  by  equation  2.9. 
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Chapter  5 


Energy  Restoring  Momentum 
Conserving  Algorithm  for  Frictionless 
Dynamic  Contact 

5.1  Introduction 

We  present  in  this  chapter  the  development  of  an  energy  restoring  algorithm  for  frictionless 
dynamic  contact.  The  objective  of  any  algorithm  is  to  represent  the  simulated  physical 
system  as  closely  as  possible.  To  ensure  this,  the  algorithmic  problem  should  inherit  the 
properties  of  the  underlying  continuum  problem. 

It  is  also  desirable  to  achieve  the  development  of  a  stable  implicit  algorithm  to  solve  friction¬ 
less  contact  problems  so  as  to  obtain  a  robust  scheme.  The  contact  algorithm  detailed  here 
has  these  conserving  properties,  in  addition  to  enforcing  the  unilateral  contact  constraint 
.and  the  velocity  constraint,  as  seen  in  Armero  &  Petqcz[1]  and  AR-MERO  ,&  PjEt6cz[2]. 
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5.2  Energy  restoring/momentum  conserving  scheme 

Since  the  goal  is  the  conservation  of  energy  and  the  momenta,  let  us  begin  by  expressing  the 
variation  of  the  total  energy  in  a  time  discretized  form  and  analyzing  the  contribution  of  the 
normal  contact  force.  The  tangent  component  which  originates  from  the  friction  phenomena 
is  not  be  included  in  the  analysis. 

5.2.1  Time  discretization  of  the  weak  form 

In  the  absence  of  external  forces  and  imposed  displacements,  but  including  the  normal  con¬ 
tribution  of  the  contact  traction,  equation  4.22  takes  the  following  form: 


en+l-£n  =  {Kn+i-}Cn)  +  {yVn+i-Wn) 

Jure 

i‘Pl+i{Yn+e^{X))  -  v>l{YnMx)))]  dr,  (5.1) 

where  we  have  left  as  an  unknown  the  instant  tn+e^  at  which  we  perform  the  closest  point 
projection  within  the  time  step  [tn,  tn+i]-  Also,  in  the  absence  of  friction,  the  contact  traction 
at  is 


(5.2) 

where  normal  at  the  contact  point  at  the  configuration  at  tn+02-  Substituting 

the  expression  of  the  contact  force  in  expressioni  5.2  into  equation  5.2,  we  obtain 


^n+l  ~  —  {^n+1  “  ^n)  +  (TVn+l  ~  ^n) 

=  f  iAr[(v’Ui(^)-¥>i(V))- 

JuFc 

ivl+iiYn+M)  -  v^2(y„^,^(x)))]  .  Un+92  dr  .  (5.3) 
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We  have  established  the  following  three  unknowns:  6i,  $2  and  We  need  to  find  expressions 
for  these  unknowns  so  as  to  satisfy  the  restrictions  imposed  by  the  conservation  properties 
we  want  the  algorithm  to  have,  namely  energy  restoration  and  linear  and  angular  momentum 
conservation. 

We  define  the  following  evolution  equation: 


9n+l  ■■=  9t  +  ^n+e2  •  [W+l(^)  -  V>l+iiYnMx)))  -  (X)))]  .  (5.4) 

In  view  of  the  above  definition,  we  define  a  new  gap  function  which  we  call  the  dynamic  gap 
function  and  denote  by  Armero  &  Petocz  [1]  [2]. 

Remark  5.1  Equation  5.4  is,  in  fact,  a  time  discretization  of  the  gap  velocity  expression, 
given  by 

^  =  1/-  [f1(x)-v2(y-(x))]  .  (5.5) 

Hence,  the  contribution  of  the  normal  contact  pressure  to  the  energy  variation  expressed  in 
equation  5.1  can  be  rewritten  as  follows: 


4+1  -  4  =  /  Mgi^,  -  gi)  dr  .  (5.6) 

Jure 

Keeping  in  mind  that  the  goal  is  the  development  of  a  stable  algorithm  in  the  form  of  an 
energy  restoring  scheme,  i.e.  one  for  which  the  energy  is  restored  to  the  system  of  bodies 
when  contact  is  no  longer  detected,  we  have  to  satisfy 


4+1  4  —  {r n+l  r n)  )  (^•'^) 

where  Vt  is,  intuitively,  the  total  regularization  potential  of  a  conservative  force  which  will 
enforce  the  unilateral  contact  constraint.  By  construction,  we  then  have 
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Pn+i  -P^  =  -  f  {U(gi^,)  -  V{gi))  dr,  (5.8) 

JuFc 

where  the  potential  U  is  the  individual  contribution  of  each  point  in  contact.  Using  equation 
5.6  together  with  equation  5.8,  we  have  the  following  expression  for  the  normal  contact 
pressure  tN  Armero  &  Petocz  [1]  [2]: 


(  U{9Ui)-U{9i) 

tjy  =  <  9n+l~9n 

1  -u-  iUsUi + s;!)) 


if  si+i  #  si  , 
if  si+i  =  si  , 


(5.9) 


where  U{g)  is  the  penalty  regularization  potential  that  is  used  to  satisfy  the  constraint 
detailed  in  equation  2.19. 

Remark  5.2  The  regularization  potential  U{g)  used  in  5.9  can  be  any  positive  definite 
function  which  is  nonzero  if  and  only  if  5^  <  0.  However,  the  function  has  to  be  such  that 
when  g^j^^  <  0  and  >  0  it  enforces  the  impenetrability  constraint  at  tn+i 


Intuitively,  we  may  imagine  the  contact  force  as  the  sum  of  the  action  of  a  series  of  “stiff’ 
springs  located  along  the  common  surface  of  contact  UTc?  where  this  force  is  a  function  of 
the  existing  gap  between  the  two  contacting  surfaces.  Armed  with  this  physical  insight,  we 
define  U (g)  as  follows: 


U{g)  :=  { 


V 


1  9 

0 


if  5  <  0  , 
if  ^  >  0  , 


(5.10) 


where  kn  is  &  (large)  penalty  parameter. 

The  rest  of  the  unknowns,  namely  9i  and  62,  are  set  so  as  to  enforce  the  conservation  of 
angular  momentum.  We  add  to  the  calculations  of  Section  3.3.1.2,  the  contribution  of  the 
normal  contact  force  to  the  variation  of  the  angular  momentum.  If  one  performs  the  standard 
substitution  in  the  variations  Stp,  one  obtains 
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[Jn+l  -  Jn]-W  =  /  tNl/n+92  '  (-X”)  X  W  -  <p^  i(Yn+9i{X))  Xw  dE 

JuPc  L  2  2  J 

=  -  V^^i(rn+ei(-X))]  XW  dr 

=  [*^71+^2  X  (v’i+i  (-x)  -  ¥5^^i(y„+0i(x)))]  -wdr .  (5.11) 


Notice  that  if  9i  =  O2  =  h 


[Jn+l  -Jn]-W  =  tN  X  •  W  dr 


=  0  Vtr  e  , 


(5.12) 


where  we  have  used  the  gap  definition  in  equation  2.17  at  and  the  unit  normal  as 
defined  by  the  closest  point  projection. 

Conservation  of  linear  momentum  is  guaranteed  by  construction.  For  completeness,  we 
include  the  proof.  Using  the  variation  and  the  results  from  Section  3.3. 1.1.,  we  obtain: 


[in+l  -^n]  ■  ®  ®] 

J'JFc 


0  Va  G  E'‘<**”*  . 


(5.13) 


Using  the  previous  results,  we  can  rewrite  5.4  as  follows  Armero  &  Pet6cz[1]  [2]: 


^n+l  =  +  ^n+i  '  [(v’«+l(^)  "  V’n+lC^n+i  W))  "  (¥>n(A-)  -  V^^(i;+i  (X)))]  (5.14) 

Remark  5.3  Equation  5.14  is  a  second  order  approximation  for  the  evolution  equation  of 
the  gap  function  in  the  continuum  setting,  expressed  in  equation  5.5. 

Remark  5.4  Notice  that  the  definition  of  the  dynamic  gap  function  gf  is  defined  by  the 
deformation  ip,  and  not  by  the  material  velocities  V,  so  that  the  restoring  scheme  can  be 
used  with  only  slight  modifications  in  a  quasi-static  problem. 
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5.2.2  Energy  restoration 

The  energy  restoration  property  of  our  scheme  for  the  frictionless  case  is  stated  in  equation 
5.7,  where  the  total  contribution  of  the  penalty  potential  is  a  positive  quantity,  that  is, 

V,  :=  [  U(gi)dr 

Jure 

>  0  .  (5.15) 

The  combination  of  the  energy  balance  expressed  in  equation  5.7  and  the  positive  definiteness 
of  the  potential  Vt  stated  in  equation  5.15  has  two  important  consequences: 

1.  For  any  time  step  tn  in  which  contact  is  detected,  i.e.  U{gn)  >  0, 

Sn<£o,  (5.16) 

where  Sq  is  the  initial  total  energy.  Hence,  the  energy  of  the  system  of  solids  will  never 
increase  during  the  time  of  contact  regardless  of  the  size  of  the  time  step  At  used  in 
the  time  marching  scheme. 

2.  For  any  time  step  t„  in  which  contact  is  no  longer  detected,  i.e.  U{g^)  =  0, 

Sn^So.  (5.17) 

Hence,  the  energy  of  the  system  of  solids  is  restored  after  release  occurs. 

We  can  conclude,  then,  that  the  proposed  scheme  is  unconditionally  (energy)  stable. 
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5.2.3  Space  discretization 

Now  we  need  to  formulate  equation  5.14  in  a  space  discretized  setting.  The  distinctive  aspect 
is  that,  when  in  contact,  we  do  not  identify  a  common  contact  surface  Ec  but  only  a  set  of 
slave  nodes  S  which  come  into  contact  with  their  master  segment.  As  explained  in  Chapter 
4  we  have  Ugiave  slave  node/master  segment  pairs.  Then,  for  each  slave  node  material  point 
Xs,  we  can  define  a  dynamic  gap  at  time  denoted  by  „  as  follows: 


9h+I  =  sin  +  -  »>S+i(?„+l(Xs)))  -  (».i(Xs)  -  «=J(f„+l(Xs)))] 

(5.18) 

With  the  notation  introduced  in  the  previous  chapter,  we  can  rewrite  equation  5.18  in  the 
following  way: 

<n+l  =  gin  +  [d;,n+:  -  l,n\  ,  (5.19) 

for  the  corresponding  displacements  of  slave  and  master  nodes  at  tn+i  and  tn-  The  initial¬ 
ization  of  the  dynamic  gap  evolution  equation  is  achieved  by  using  the  (real)  gap  function 
gs,n  (be.  the  one  obtained  by  performing  a  closest  point  projection  procedure  at  t„,  the 
last  time-step  before  contact  is  detected).  Equation  5.18  is  a  second  order  approximation  of 
equation  2.42  for  the  evolution  of  the  real  gap  gs  in  a  space  discretized  formulation. 

Remark  5.5  Note  that,  in  the  one  dimensional  case,  the  dynamic  gap  g^  and  the  (real)  gap 
gs  coincide,  as  the  approximation  is  given  by  the  change  in  the  normal  Ug. 

The  contribution  of  the  normal  contact  force  to  the  nodal  contact  vector  is  then  given  by 


(5.20) 


where  is  defined  in  equation  5.9. 
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5.2.4  Contact /release  conditions 

An  independent  part  of  a  contact  algorithm,  in  addition  to  a  contact  detection  algorithm 
and  the  actual  contact  force,  is  the  appropiate  choice  of  a  contact /release  criterion.  The 
contact/release  logic  developed  for  this  scheme  is  presented  in  Table  5.1,  below. 

As  mentioned  previously,  the  dynamic  gap  „  is  initialized  with  the  value  of  the  “real” 
gap  from  the  converged  values  of  the  previous  time  step  before  contact  is  detected  at  tn+i- 
Contact  is  detected  at  if  is  negative,  which  is  our  indication  that  penetration  has 
occurred.  Subsequently,  contact  is  detected  if  the  dynamic  gap  gg^n+i  is  negative.  As  the 
contact  pressure  tjv  depends  on  the  contact  states  at  tn+i  and  tn,  it  only  vanishes  when  both 
of  these  states  are  released  states.  When  contact  is  no  longer  detected  at  t„+i,  we  observe 
that  tjv  becomes  negative;  this  represents  the  final  “kick”  that  restores  the  energy  back  to 
the  system  of  bodies. 

An  important  aspect  of  our  algorithm  is  that  we  enforce  the  unilateral  constraint  at  tn+i- 
At  the  first  time  step  when  contact  is  detected,  U{gg  „)  =  0,  so  the  penalty  regularization  is 
applied  to  enforce  g^+i  =  0-  Many  standard  contact  algorithms  impose  the  gap  constraint 
at  tn+§)  which  leads  to  an  oscillatory  behaviour  of  the  gap  function. 

5.3  Enforcement  of  the  gap  velocity  constraint 

Different  mechanical  contact  problems  involve  different  contact  times.  By  contact  times  we 
mean  to  the  period  of  time  two  bodies  may  remain  in  contact.  This  contact  time  depends 
on  the  mechanical  properties  of  the  material  and  the  relative  velocity  between  the  bodies 
before  collision. 

In  some  cases,  we  may  want  to  investigate  the  contact  problem  with  very  small  time  steps 
as  we  might  be  interested  in  what  happens  during  the  contact  time.  In  other  cases,  we 
may  be  more  interested  in  the  collision  paths  and  rebound  velocities  so  we  may  perform  the 
simulation  with  bigger  time  steps. 
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_ Table  5.1:  Contact/release  logic. _ 

Let  cents, „  =  contact  flag  at  tn  (.true,  or  .false.),  and  gs,n+i  the 
(real)  gap  at  tn+i  for  slave  node  S.  Then,  conts,n+i  is  defined  as  follows; 

IF  (conts,„  .or.  -le-O))  THEN 

Compute  gf^n+i  using  5.18 
IF  «n+i-le.0)  THEN 
conts,„+i  =  .true. 

ELSE 

cents, „+i  =  .false. 

IF  (^s,n+i.ge.0)  THEN 

The  dynamic  gap  will  be  initialized  to  the  curren 
9s,n+i  when  evaluating  5.18  in  the  next  time  step 
ELSE 

The  dynamic  gap  will  be  initialized  to  the  current 
9s,n+i  when  evaluating  5.18  in  the  next  time  step. 

ENDIF 

ENDIF 

ELSE 

conts,n+i  =  .false. 

ENDIF 
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The  algorithm  presented  works  for  any  case  including  the  two  extreme  cases  detailed  above. 
If  the  contact  time  is  “big”  compared  to  the  time  step  chosen  to  perform  the  simulation,  then 
the  simuation  falls  into  the  category  of  persistent  contact.  As  discussed  above,  during  persis¬ 
tent  contact  we  encounter  the  gap  velocity  constraint  as  stated  in  equation  2.31.  Satisfaction 
of  this  constraint  has  been  researched  by  Lee  [21]  and  by  Taylor  &  Papadopoulos  [33]. 

The  goal  in  this  section  is  to  modify  the  previously  presented  algorithm  to  satisfy  this 
additional  constraint,  while  maintaining  the  conservation  properties  achieved  so  far.  To  this 
end,  we  restate  the  weak  form  of  the  balance  of  momentum  equation  presented  in  equation 
4.3,  in  the  following  form: 


E  +  ^  :  Grad((5^«)  dQ 

=  f  [tNSg  -  tTpS^^]  dP.  (5.21) 
Jure 

Also,  equation  4.8  can  be  rewritten  as  follows: 


for  bodies  i  =  1, 2. 


At 


(5.22) 


We  use  the  following  approximation  to  the  gap  velocity  expression  given  in  equation  2.27: 

ht  =  [v/(x)  -  (i7(x))]  •  ut ,  (5.23) 

for  any  t  G  [t„,  tn+i],  where  we  have  used  the  notation 


^  ^(Ln)  (5.24) 

We  perform  a  penalty  regularization  on  the  linear  momentum  p  for  the  points  belonging 
to  the  neighbourhood  of  the  common  boundary  0{rc)  to  enforce  the  constraint  on  the  gap 
velocity  g.  Hence,  we  present  the  definition  of  the  regularized  linear  momenta  as  follows: 


52 


5.3.  ENFORCEMENT  OF  THE  GAP  VELOCITY  CONSTRAINT 


pI  =  +  VxeO(rc),  (5.25) 

=  p^v^  -  Pp,thti^t  Vy(x)eO(rc),  (5.26) 

where  p\  and  denote  the  modified  linear  momenta  for  each  material  point  X  in  the 
neighbourhood  of  the  surface  of  body  1  and  its  corresponding  projection  y  (x)  onto  body  2, 
respectively. 

Remark  5.6  This  modification  of  the  definition  of  the  linear  momenta  in  the  neighbourhood 
of  Fc  leads,  in  the  finite  element  implementation,  to  the  use  of  a  lumped  mass  matrix  instead 
of  a  consistent  mass  matrix. 

The  penalty  Pp^t  is  defined  in  the  following  way: 


where  rrip  is  a  (large)  penalty  parameter. 


if  U{gt)  #  0  , 

otherwise  , 


(5.27) 


The  modification  to  the  definition  of  the  linear  momenta  detailed  in  equations  5.25  and 
5.26  and  the  subsequent  change  in  the  dynamic  formula  5.22  remains  active  only  while  the 
impenetrabilty  constraint  {gf  =  0)  is  active,  as  expressed  in  equation  5.27.  Notice  that 
for  interior  points,  the  definition  of  the  linear  momenta  remains  unchanged  and  the  time 
updating  formula  4.8  remains  valid.  This  modification  enables  us  to  introduce  a  penalization 
on  the  gap  velocity  g  in  an  algorithmic  way. 


This  enhancement  to  the  scheme  to  enforce  the  gap  velocity  constraint  preserves  both  its 
conservation  and  restoring  properties.  We  prove  this  statement  in  the  following  section. 


5.3.1  Properties  of  the  proposed  scheme 

Our  goal  has  been  to  enforce  the  constraint  on  the  gap  velocity,  as  expressed  in  equation  2.31, 
while  maintaining  the  validity  of  the  conserving  properties  of  our  scheme.  Next,  we  inves- 
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tigate  these  properties  and  prove  that  we  have  created  an  energy  restoring  and  momentum 
conserving  scheme  with  the  added  characteristic  of  satisfying  the  gap  velocity  constraint. 
Proposition  5.1 

In  the  absence  of  external  forces  and  imposed  displacements,  we  have  for  a  typical  time 
interval 


1.  The  linear  momentum  is  conserved,  that  is: 

Ln+l  —  I'n  j 

2.  The  angular  momentum  is  conserved,  that  is: 

•In+l  ~  Jn  ) 


3.  The  energy  evolves  as 


En+l  +  Pn+l  +  A4n+1  =  £n  +  Pn  A  A4n  , 


where 


Mt  := 


> 


dQ 


(5.28) 


(5.29) 


(5.30) 


(5,31) 


Proof: 

1.  We  prove  that  the  modification  to  the  linear  momenta  does  not  contribute  to  the  total 
linear  momentum  of  the  system.  Therefore,  the  conservation  of  linear  momentum 
proved  in  the  previous  section  still  holds. 


u  = 


0. 


(5.32) 
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2.  The  conservation  of  the  angular  momentum  is  proved  in  a  similar  way.  We  first  prove 
that  the  contribution  to  the  total  angular  momentum  of  the  system,  by  the  modification 
to  the  linear  momenta,  vanishes.  That  is. 


-u 


f  ^tX  P°‘Vt  df2+  f  -  v>l{Yt{X))]  X  Pp^thti't  dQ 

Jo{rc) 


X  P°^Vt  dQ  +  /  Qtvt  X  ppj^htut  dQ 
i  Jni')  Jo(r^) 


=  Jt, 


(5.33) 


where  the  integrand  of  the  surface  integral  vanishes  due  to  co-linearity. 

As  before,  the  conservation  of  this  new  quantity  still  holds,  so  the  conservation  of  the 
total  angular  momentum  is  confirmed. 


3.  In  terms  of  the  energy,  we  know  that  the  following  property  is  preserved: 


^n+l  +  yVn+1  +  'Pn+l  =  +  VVn  -I"  Vn  ,  (5.34) 

where  we  denote  by  ICt  the  modified  total  kinetic  energy,  stated  as 

■=  ^  ^  p?  dQ  ,  (5.35) 

for  t  €  \tn,  tn+i]-  After  some  algebraic  manipulation,  we  arrive  at  the  following  expres¬ 
sion: 


where 


K,= 

lCt  + 

/  Pp,tht 

i+Mj 

M  1 Y 

JoiFc) 

2  ' 

KP^  P^J. 

Ml 

:= 

1  Pp,t^t 
Jo{rc) 

X  1 

2 

> 

0 

dQ  , 


(5.36) 


dQ 


(5.37) 
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5. 3. 1.1  Energy  restoration 

The  energy  restoring  property  still  holds  due  to  equation  5.36  and  the  fact  that  the  added 
kinetic  energy  A4t  introduced  by  the  penalty  regularization  on  the  linear  momentum  is 
always  positive  and  vanishes  at  a  released  state. 

5.3.2  Finite  element  implementation 

The  previously  stated  modification  to  the  definition  of  the  linear  momentum  of  the  problem 
is  easily  implemented  within  a  finite  element  formulation.  To  this  end,  we  rewrite  the  linear 
momentum  for  a  typical  slave  node/master  segment  as  follows: 


Ps,t  —  “1“ 

Ms,L^s,t  d"  )  (5.38) 

where  t  €  a  typical  time  increment,  and  we  denote  by  hs,t  the  discrete  counterpart 

of  the  gap  velocity  definition  2.42.  It  is  defined  as 


^‘'master 

1=1 


MI 

t 


(5.39) 


We  have  used  in  equation  5.38  the  lumped  mass  matrix  Ms,l  for  the  contact  element,  i.e. 
the  slave  node  and  the  master  nodes  that  participate  in  the  master  segment.  It  is  given  by 


:= 


“^dim 


G  ^  (5.40) 
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to  simplify  the  notation.  In  equation  5.38,  nis^t  denotes  a  penalty  mass  added  to  the  con¬ 
tacting  slave  and  master  nodes,  which  depends  on  the  contact  states,  as  follows: 


m 


Syt 


1 


mp 

0 


if  tNs,t  ^  d  , 
otherwise  , 


(5.41) 


where  nip  is  a  large  penalty  parameter.  As  nip  oo,  the  constraint  hg^n+i  =  0  for  a  typical 
interval  [tn,tn+i]  of  contact  is  enforced. 


Rewriting  equations  4.42  and  4.43,  and  substituting  the  linear  momentum  p  for  our  new 
intermediate  variable  p  defined  in  equation  5.38,  we  have  the  following: 


(ds,n+l  "h  r)^s,L  ip^s,n+lhs,n+l^s,n+l  F  '^s,nhs,n^s,n}  ?(b.42) 


?(”+5)  r  s(^+5) 


?("+2) 


\Vs,n-^l  ^5,nj  fs,mt  f s,{c, mass)  ^  ext  » 


(5.43) 


where  we  denote  by  the  contact  force  given  by 


and 


-{n+|)  _  ;(n+5)  ^  -5(n+|) 

I s,(c,mass)  **  > 


(5.44) 


*  • —  A,  {p^s,n+lhs,n+l^$,n+l  'l^s,nhs,n^s,n^ 


At 


(5.45) 


This  added  force  enforces  the  constraint  on  the  gap  velocity,  and  vanishes  for  the  nodes  that 
are  not  in  contact. 


5.4  Addition  of  energy  dissipation 

Impact  phenomena  are  difficult  to  simulate  due  to  the  high  frequency  modes  present  in 
the  motion  of  the  body  during  and  after  the  impact.  The  need  to  introduce  numerical 
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dissipation  into  the  algorithmic  scheme  becomes  apparent  during  short-term  simulations 
where  the  interval  in  which  the  bodies  are  in  contact  is  long  compared  to  the  time  step  At. 

To  this  end,  we  introduce  a  modification  to  the  previous  energy  restoring  algorithm.  During 
persistent  contact  for  the  time  step  [tn,tn+i],  )  ^  0  and  U{gU  ^  0.  Thus,  the 

normal  contact  pressure  tff  takes  the  following  form  for  the  potential  function  described  in 
equation  5.10: 


tw 


{qU+iY  -  (gUf 

- TTd - 

-l>^N  {gin+l  +  gin)  ■ 


We  perform  a  modification  of  the  previous  expression,  as  follows: 


(5.46) 


tN  =  -Kn  W,«+1  +  (1  -  ’’)«?,») 


(6.47) 


for  >  5. 


Remark  5.7  At  the  first  and  last  time  steps,  the  expression  for  the  contact  pressure  remains 
unchanged  with  respect  to  the  energy  restoring  one,  so  that  we  satisfy  the  impenetrability 
constraint  at 


We  now  prove  that  we  have  energy  dissipation  for  any  >  |.  Let  us  consider  the  time 
interval  [t„,  where  we  have  persistent  contact  (i.e.  contact  is  detected  both  at  tn  and 
at  tn+i).  We  have  the  following  expression  for  the  variation  of  the  total  energy  5^: 
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^n+l  Ffi  —  tpf  (^n+l  “■  fi'n) 

JPc 

=  -  [  kn  i'&gt+i  +  (1  -  d)gi)  (5^+1  -  g^)  dF 

J  Pc 

=  5  W+i  +  9«)  +  (>’ -  5)  (Ai  -  (sLi-si)dr 

=  -  /  (uigLi)  -  u{9i))  dr 

J  Pc 

~  “  5) 

=  -  ('Pwi  -  P«)  “  “  5)  “  Sif  dr  .  (5.48) 

Remark  5.8  Notice  that  we  revert  to  the  energy  restoring  scheme  when 'd  =  ^. 


59 


Chapter  6 


Frictional  Dynamic  Contact 


6.1  Introduction 

A  relevant  feature  of  the  friction  phenomenon  is  the  presence  of  energy  dissipation,  in  the 
form  of  heat.  Though  from  the  algorithmic  point  of  view,  this  dissipation  is  advantageous 
(provided  the  algorithm  inherits  this  dissipative  property),  the  treatment  of  friction  becomes 
complicated  due  to  the  discontinuous  slip/stick  behaviour  of  Coulomb  friction. 

Following  the  characteristic  line  of  our  research,  our  aim  in  this  part  of  the  work  is  to 
design  a  friction  algorithm  that  will  inherit  the  physical  properties  of  the  underlying  contin¬ 
uum  problem.  A  dissipative  algorithm  not  only  simulates  the  actual  physical  process,  but 
also  guarantees  stability  of  the  algorithm  in  the  sense  of  energy  control  (see  Armero  & 
Pet6cz[4]). 

6.2  Dissipative  friction  algorithm 

As  stated  above,  friction  is  a  dissipative  phenomenon  so  it  is  of  interest  to  study  the  evolution 
of  the  energy  of  the  system  in  the  presence  of  frictional  contact  forces. 

To  investigate  the  contribution  of  the  tangential  contact  forces  to  the  variation  of  the  total 
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energy  £t,  we  state  the  variation  of  the  energy  for  a  typical  interval  of  time  [tn,tn+i]  as  in 
equation  5.1,  where  we  have  used  (pn+i  —  as  our  admissible  variation  5(p'm.  the  weak  form 
equation  3.1.  By  performing  this  substitution  in  the  expression  for  5^^  in  equation  3.4,  the 
energy  variation  now  takes  the  form 


■'n+1 


-£n 


-L 

-I 


tNigUi  -  oi)  dr 


^Ta^a(5 


ll^’n+l  -  ¥>nl 


(n+l 


+  ^n+i  ( 


^n+1,/3 


*^71+1 


diC6-l) 


Analogous  to  the  definition  of  the  dynamic  gap  function  g‘^,  let  us  define  the  dynamic  slip, 
denoted  by  in  the  following  way  [3]; 


(n+l) 


“1”  {,^n+l,0  ^n,^) 


*^n+|  > 


(6.2) 


71^  JL 

where  the  geometrical  quantity  A^^^  is  defined  in  equation  2.44  and  is  evaluated  at  the 
configuration  at  Also,  the  definition  of  the  normal  is  obtained  by  means  of  a 

closest  point  projection  performed  at 


Remark  6.1  Equation  6.2  is  a  second  order  approximation  of  the  evolution  equation  for 
the  tangential  slip  expressed  in  equation  2.43. 

Substituting  definition  6.2  into  equation  6.1,  the  expression  for  the  energy  variation  becomes 


^n+l  £n  —  j ^  5n)  ^Ta  ^^n+1  dF  .  (6-3) 

Remeirk  6.2  Notice  that  the  definition  of  the  dynamic  slip  in  equation  6.2  involves  the 
deformation  (p  and  not  the  material  velocities  V,  so  that  this  contact  quantity  can  be  used, 
with  only  slight  modifications,  in  a  quasi-static  problem. 
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6.2.1  Split  operator  integration  algorithm 

Recall  the  equations  that  govern  the  friction  phenomenon,  namely, 


<t> 

•=  ll*r||  —  <  0  , 

(6.4) 

7 

>  0, 

(6.5) 

7(/. 

=  0, 

(6.6) 

7<^ 

=  0, 

(6.7) 

and  the  regularized  slip  rate  relation,  given  by 

On  the  other  hand,  contact  kinematics  yields  the  following  expression  for  the  tangential  slip 
rate  vt' 


VT  =  .  (6.9) 

We  substitute  the  dynamic  slip  rate  for  in  equation  6.9  [3],  to  obtain  the  following 
expression: 


VT  =  .  (6.10) 

Next,  the  slip  relation  6.8  combined  with  equation  6.10  is  discretized  in  time  by  means  of  a 
generalized  mid-point  approximation  of  the  form 


-  d'")  =  ^ 


(6.11) 
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for  any  0  €  (0, 1].  We  have  used  an  approximation  to  the  reference  metric,  denoted  by  Map', 
the  nature  of  this  approximation  is  explained  below. 

The  slip  surface  function  (p  and  Kuhn- Tucker  conditions  for  a  typical  time  interval 
now  take  the  following  form: 


0  = 

ll^rn+sl  -  fJ-tN  <  0  , 

(6.12) 

l> 

IV 

0, 

(6.13) 

A7  (j)  = 

0. 

(6.14) 

To  integrate  equation  6.11,  while  enforcing  the  Kuhn-Tucker  complementary  conditions  6.13 
and  6.14,  a  operator  split  strategy  is  used: 

1.  Trial  state 

As  in  the  standard  return  mapping  algorithm,  we  first  evaluate  the  trial  state  by 
assuming  a  stick  phase  by  setting  A7  =  0,  thus  obtaining  an  expression  for  the  trial 
state  tangent  stress  at  tn+9,  that  is, 

=  ,  (6.15) 

where  we  have  assumed  that  stick  conditions  apply  for  N  time  increments.  In  other 
words,  A7  =  0  for  n  =  0,  •  •  •  ,  W  —  1  and  the  initial  point  of  contact  at  the  initiation 
of  the  N  stick  time  steps  is  f  =  |^. 

To  decide  which  of  the  frictional  states,  slip  or  stick,  is  present,  we  evaluate  the  function 
^  using  this  trial  stress.  We  denote  this  value 

2.  Stick  phase  ■■■  -  ■  ■ 

If  <  0,  then  A7  =  0  must  hold  to  satisfy  the  condition  6.14.  Consequently,  we 
satisfy  the  stick  condition  and  the  trial  tangent  stress  in  equation  6.15  is  the  actual 
frictional  stress,  i.e.. 
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4.  _  j.zriat 

^Tn+e  —  ^Tn+$ 


(6.16) 


3.  Slip  phase 

Assume  that  step  iV  +  1  is  a  slip  step,  that  is 


>  0  ,  (6.17) 

meaning  that  we  need  to  correct  the  trial  frictional  stress  so  that  it  satisfies 

condition  6.14.  Since  a  slip  state  implies  A7  >  0,  we  know  that 

<^  =  0  .  (6.18) 

Therefore,  substituting  the  expression  of  the  trial  state  frictional  stress  in  equation 
6.15  into  equation  6.14  yields 


w.  -  w  [1  +  j  , 

implying  that  the  trial  and  actual  stress  states  are  parallel,  i.e. 


(6.19) 


Arial 

^TN+e 


Arial 

^TN-^e 


txN+e 

II^TAT+^II 


(6.20) 


This  last  result  is  crucial  to  solve  the  problem  in  closed  form.  Using  equation  6.20 
together  with  condition  6.18,  we  obtain 


tTN+9  =  IJ'tN 


Arial 

^TN-^e 


lt 


trial  5 
TN-\-e 


(6.21) 


withA7=i£;. 


A  slip  frictional  state  involves  relative  slip  between  the  two  surfaces  and  we  need  to 
update  our  stick  point  to  the  new  contact  point.  That  is, 
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Remark  6.3  The  proposed  frictional  scheme  uses  a  constant  reference  metric  Map  evaluated 
at  the  stick  point  for  simplicity.  As  kt  >  0,  the  slip  is  very  small  during  a  stick  step,  so 
this  simplification  does  not  add  any  significant  limitation  to  the  frictional  scheme. 


6.3  Energy  evolution  of  the  frictional  algorithm 

In  this  section  we  prove  that  our  scheme  is  unconditionally  dissipative.  Including  the  results 
from  the  previous  chapter  for  the  energy  evolution  in  the  presence  of  normal  contact  forces, 
we  have 

(^n+l  +  'Pn+l)  -  i^n  +  'Pn)  =  —  J  txa  (^n+1  ■“  •  (6.22) 

For  each  frictional  state  we  have  a  different  expression  for  the  tangent  friction  force.  We 
treat  each  one  individually: 

1.  Stick  phase 

As  in  the  previous  section,  assume  that  there  are  N  consecutive  steps  which  satisfy 
the  stick  condition.  For  later  use,  we  define  the  following  notation: 

Vt  ■■=  - 1''  ■  (6.23) 

Also,  we  define  a  norm  by  means  of  the  approximate  metric  Map: 


WvWm  =  rfMapn^  .  (6.24) 

Using  expression  6.16,  the  change  in  total  energy  for  a  typical  time  step  [t„,  tn+i]  is  as 
follows: 
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(Sn+l  +  Vn+l)  -  {Sn  +  Vn)  =  -  Kt  (Cnfi  -  Ma0  dT 

=  M^0 

(^(^nfi  -  -  (1  -  e)&  -  dr 

=  -  [  I^Avi+l  -  Vn'')Mai3 
Jure 


=  -/ 

JuFc  L  \ 


h[ 


d  2 
>1+1  M 


1. 


-M^-2)Nn+l-^X 


dr . 


(6.25) 


The  variation  of  the  total  energy  after  N  steps  is  then 


x;  (£„+■  +  P„+i)  -  (£.  +  P„)  =  -  f  dr 

^  n^oduFc 

<  0  ,  (6.26) 

where  t}^  =  0,  since  the  stick  point  was  initialized  as  ^q-  This  observation,  together 
with  the  fact  that  9  >  allows  us  to  establish  that  the  two  integrands  in  6.26  are 
positive. 

Remark  6.4  It  is  important  to  notice  that  even  though  the  scheme  introduces  some 
dissipation  in  the  stick  phase  (due  to  the  approximate  way  in  which  we  satisfy  the  stick 
constraint),  by  choosing  0  =  |,  we  can  minimize  this  artificial  dissipation.  Also,  notice 
that  in  the  limit  as  /cy  -»  oo,  this  dissipation  vanishes. 

2.  Slip  phase 

To  calculate  the  variation  of  the  energy  for  a  slip  step,  we  use  expression  6.21  to 
substitute  into  equation  6.22,  and  we  obtain 
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{Fn+l  +  Vn+l)  —  {En  +  Fn)  —  “ 


^trial 

tii  T(xn-\-6 
atrial  |j 

rTn+eW 


{C,  ~  ii") 


_  f  f  ,cd,a  _  td,a\  _  (cd,a  _ 

Mo^  {Hdii  -  -  {1  - - f-'*))  dr 


Wttrial  II  Wn+1  Vn 


UFc  11*^71+^11 

2 


+  {o-^)iviii  +  vi’^)]  dr 


f  r  I^TfitN  fhi+i\?M  ll’?nllA/\ 

UmiiewK  2  2  j 

-S^A^-\)hU-n‘A  dr.  (6.27) 

IrTn+ell  ^  J 

If  we  now  consider  the  contributions  of  the  N  stick  time  steps  followed  by  a  final  slip  time 
step,  we  obtain  the  following  expression  for  the  variation  of  the  energy  for  the  time  interval 

[0,  tjv]: 


(5jv  +  Vn)  ~  {Fq  +  Vo 


'!  <_/'  hff+iWli 

/  —  /  ±tTial  II  9 

«/urc  L  ^Tn+^ll  ^ 


—  Kt 


l|2-  1 

2 


dr 

(6.28) 


where  we  have  discarded  the  ||77n+i  —  ^jfiH  terms  to  arrive  at  the  inequality.  The  first  term  of 
the  integrand  is  positive  due  to  the  fact  that  normal  contact  involves  compressive  tractions 
on  the  contact  surface,  that  is,  >  0;  and  the  expression  in  parenthesis  of  the  second  term 
is  positive  due  to  the  inequality  6.17  which  holds  for  the  slip  time  step. 

Thus,  the  energy  is  dissipated  throughout  frictional  contact  conditions;  this  contributes  to  the 
numerical  stability  of  the  algorithm.  It  can  be  shown  that  the  conventional  frictional  schemes 
cannot  be  proven  to  dissipate  unconditionally  if  the  dynamical  slip  is  not  introduced  into 
the  scheme;  this  not  only  makes  the  algorithm  energetically  unstable  but  also  represents  an 
infeasible  physical  situation. 
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6.4  Space  discretization 

We  now  develop  the  expressions  for  the  frictional  contact  problem,  in  a  finite  element  for¬ 
mulation  in  the  two  dimensional  setting. 

The  tangent  traction  is  expressed  in  a  space  discretized  setting.  For  each  slave  node  S  with 
coordinates  Xs  we  can  define  the  variation  in  the  following  way: 


pdjQ:  _ I 


mil 


n+j 


9s,n+^ 
- ; — -D 


n+i 


(dn+1  5 


=:frn+i 


with  the  vectors  and  defined  as: 


0 

1 

+ 

_ 1 

and  T_  ,  1  = 

^+2 

(6.29) 


We  have  used  the  following  compact  notation  to  express  the  displacements  of  the  slave  and 
master  segment  nodes  {S,  Ml,  M2, ....}: 


d\i 

d^\ 

n+2 

jAf2 

n+f 


(6.30) 


We  assume  that  friction  is  present  in  the  problem  as  long  as  the  bodies  are  in  contact, 
so  friction  will  be  added  to  the  problem  following  the  same  contact/release  criteria  of  the 
frictionless  problem  (see  Chapter  5). 

Using  the  previous  notation,  the  total  contact  force  in  equation  5.20  has  the  following 

form: 
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fs,c 


(6.31) 


Following  standard  procedure  in  the  finite  element  formulation,  the  force  contribution  of  each 
contact  element  is  assembled  as  in  equation  4.36.  The  expression  for  the  frictional  traction 
It  for  a  particular  contact  element  {S,  Ml,  M2, . . . }  for  a  typical  time  interval  [tn,  tn+i]  takes 
the  following  form: 


where 


tTn+Oa  — 


ftrial 

if  >  0 


(6.32) 


+  (1  - 


(6.33) 


The  stick  point  takes  its  value  from  the  converged  value  of  the  slip  from  the  last  time  step 
for  which  the  slip  condition  was  detected  (i.e.  for  which  ftrial  >  0). 
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Chapter  7 


Sorting  Algorithm 


7.1  Introduction 

One  of  the  most  important  parts  of  a  good  contact/impact  algorithm  is  the  contact  detection 
scheme,  in  terms  of  its  speed  and  efficiency.  The  contact  detection  becomes  a  cumbersome 
task  as  the  system  of  bodies  becomes  large.  This  task  would  include  performing  closest 
point  projections  among  N  bodies.  So  the  need  arises  for  a  more  efficient  way  of  discover¬ 
ing  the  geometric  relationships  among  objects  within  a  working  space,  because  this  process 
accounts  for  a  significant  portion  of  the  computational  effort  of  solving  contact/impact  prob¬ 
lems  among  many  bodies. 

Contact  detection  can  be  defined  as  finding  the  members  of  a  set  of  points  that  lie  inside 
a  subregion  of  an  N  dimensional  space.  In  our  formulation,  the  set  of  points  is  the  set  of 
nodes  which  lie  on  the  surface  of  each  body,  and  each  body  defines  a  subregion  in  a  two 
dimensional  space. 

In  a  multi-body,  system,  it  makes  sense  to  distinguish  two  phases  in  the  contact  detection 
procedure:  a  spatial  sorting  phase  and  a  contact  resolution  phase  (see  WILLIAMS  &  O’- 
Connor[35].  The  spatial  sorting  enables  us  to  find  all  the  pairs  of  bodies  which  could  be 
potential  contactors  and  the  contact  resolution  phase  finds  the  actual  two  points  on  the 
surface  of  each  pair  of  bodies. 
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In  this  chapter  we  expand  on  the  issues  involved  in  the  spatial  sorting  phase  of  the  contact 
detection  algorithm. 


7.2  Overview  of  various  sorting  techniques 

In  general,  it  is  necessary  to  consider  a  system  of  bodies  which  evolves  with  time.  The 
efficiency  of  the  contact  detection  algorithm  can  be  greatly  improved  if  one  can  assume  a 
priori  knowledge  of  how  this  system  will  evolve. 

Algorithms  assuming  a  priori  knowledge  can  be  extremely  efficient  but  have  the  disadvantage 
that  the  range  of  problems  they  can  tackle  is  limited.  The  following  is  a  brief  overview  of 
the  kind  of  a  priori  knowledge  that  can  be  used  to  produce  a  fast  detection  algorithm  taken 
from  Williams  &  0’Connor[35]: 

1.  Fixed  topology 

Fixed  topology  can  be  found  in  finite  element  algorithms  when  the  relative  position  of 
the  elements  remains  unchanged  during  the  simulation  of  a  particular  problem. 

2.  Slowly  varying  topology 

The  objects  move  around  only  a  small  amount  according  to  some  metric,  so  that  each 
object  only  interacts  with  its  neighbours  and  we  only  keep  track  of  a  small  number  of 
possible  contactors  per  body.  If  one  can  keep  track  of  the  characteristic  velocities  in 
a  problem,  then  it  is  also  possible  to  check  for  contact  after  only  a  certain  number  of 
time  steps  instead  of  after  every  one. 

3.  Spatially  sparse  systems 

If  the  system  is  very  sparse,  it  makes  sense  to  project  trajectories  in  such  a  way  that  we 
only  check  for  the  intersection  of  cones  or  cylinders  in  a  space-time  system.  Animation 
applications  have  used  this  scheme  successfully,  (see  e.g.  Dworkin[9]). 

4.  Exhaustive  spatial  sorting  schemes 

If  the  scheme  makes  no  a  priori  assumptions  about  the  evolution  of  the  problem  and 
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it  is  based  only  on  the  present  geometric  configuration,  it  is  called  exhaustive.  These 
schemes  are  more  general  and  robust  but  are  potentially  slower  than  non-exhaustive 
schemes. 

5.  Spatial  sorting  algorithms 

Spatial  sorting  gives  us  a  valuable  tool  to  decide  which  bodies  should  be  considered 
for  a  more  detailed  contact  resolution.  It  seeks  to  avoid  the  all-to-all  body  search  for 
contact  at  each  time  step.  For  a  small  quantity  of  objects,  this  all-to-all  method  is 
acceptable.  (Shi  considers  a  maximum  of  64,  in  [28].)  But,  for  example,  a  simulation 
involving  exhaustive  searching  over  a  thousand  objects  can  become  computationally 
prohibitive. 

In  this  work  we  use  a  spatial  sorting  technique  combined  with  the  traditional  closest  point 
projection  as  our  contact  resolution  method.  Below,  we  review  several  of  the  most  powerful 
methods  as  classified  by  Williams  &  0’Connor[35]  and  in  subsequent  sections  we  describe 
in  some  detail  the  binary  space  partitioning  which  has  been  used  in  the  present  work. 

1.  Grid  subdivision 

The  grid  subdivision  method  discretizes  in  a  uniform  way  the  simulation  volume  or  area 
into  rectilinear  cells,  each  cell  enclosing  one  or  more  objects.  Objects  are  associated 
with  each  cell  based  on  their  coordinates  and  neighbouring  objects  are  detected  by  their 
cell  assignment.  The  performance  of  this  method  depends  greatly  on  the  homogeneity 
of  the  spatial  distribution  of  the  objects  within  the  working  space  and  is  not  a  good 
overall  methodology  for  a  wide  range  of  problems. 

2.  Adaptive  grid  methods 

Adaptive  cell  methods  are  used  to  avoid  the  problems  associated  with  simple  grid  sub¬ 
division,  though  at  the  cost  of  managing  multiple  cell  dimensions.  With  this  approach, 
the  scheme  suffers  when  the  distribution  of  objects  becomes  homogeneous. 

3.  Body  based  cells 

Another  alternative  is  to  base  the  cell  surrounding  an  object  on  its  centroid.  Any 
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object  lying  within  this  cell  is  considered  to  be  a  potential  contactor  and  added  to  a 
neighbour  list.  The  multi-pole  algorithm  of  Greengard  &  Rohklin  [11]  makes  use 
of  a  similar  approach  in  the  area  of  molecular  dynamics. 

4.  Spatial  heap-sort 

Heap-sort  is  one  of  several  algorithms  used  to  sort  the  objects  into  an  ordered  list 
which  at  the  same  time  gives  an  indication  of  the  location  of  each  particular  object. 
In  this  case,  the  key  to  the  ordered  list  is  the  object’s  coordinate  along  one  or  more 
global  axes.  This  method  has  been  applied  to  particle  hydrodynamics  by  Swegle 
[32].  A  modified  heap-sort  algorithm  called  DFR  has  been  developed  by  Williams  & 
O’Connor  [35]  and  applied  to  baseline  granular  simulation  problems. 

5.  Tree  methods 

Binary  sort/search  algorithms  provide  a  flexible  and  general  methodology  for  contact 
detection  in  two  dimensional  problems.  The  octree  sort /search  algorithm  was  derived 
from  the  binary  one  to  handle  problems  in  three  dimensions.  The  technique  consists  of 
treating  the  objects  as  rectilinear  cells;  however  in  this  case  the  cells  are  ordered  into 
a  tree  data  structure.  To  ensure  the  creation  of  a  balanced  tree,  Knuth  [17]  suggests 
that  insertion  into  the  tree  should  be  done  in  a  relatively  random  fashion.  The  time 
required  to  create  a  binary  tree  is  of  the  order  0{Nlog{N)),  where  N  is  the  total 
number  of  objects  in  the  problem;  idem  for  the  traversing  of  the  tree.  In  short,  such 
an  algorithm  is  a  valuable  improvement  over  an  all-to-allsearch,  which  is  0{N^). 

In  our  work,  we  implement  a  binary  tree  database  structure  to  perform  the  sorting  (see  e.g. 
Munjiza  et  AL[24]  and  Bonet  &  Peraire[6].  It  has  proven  to  be  a  versatile  tool  that 
can  be  applied  to  a  wide  range  of  problems  in  the  area  of  contact/impact  simulations. 
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level  0 


level  1 


level  2 


level  3 


Figure  7.1:  Schematic  drawing  of  a  binary  tree 


7.3  Binary  space  partitioning 

7.3.1  Binary  tree  structure 


Originally,  tree  structures  were  developed  to  store  in  a  systematic  way  a  collection  of  data  in 
order  to  enable  quick  access.  A  tree  consists  of  nodes  where  the  actual  data  is  stored.  Each 
node  is  extended  by  the  addition  of  two  links  to  two  other  nodes  known  as  the  left  child  and 
the  right  child.  The  node  from  which  a  particular  node  springs  is  called  the  parent  node. 
Each  tree  has  a  starting  node  which  we  name  root.  Also,  at  each  node,  there  is  a  subtree 
originating  from  it  so  that  the  node  becomes  the  root  of  this  subtree. 

This  definition  establishes  a  hierarchy  of  nodes:  the  root  at  the  top  level;  0,  1  or  2  nodes  at 
the  next  level,  each  of  which  in  turn  has  0,  1  or  2  nodes  at  the  next  level  of  hierarchy;  and 
so  forth.  A  node  without  children  is  said  to  be  a  leaf.  This  hierarchical  structure  inspires 
the  graphical  representation  shown  in  Figure  7.1. 
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7.3.2  Space  partitioning  using  a  binary  tree 

The  first  step  in  a  sorting  algorithm  is  to  construct  the  binary  tree.  Consider  a  problem 
with  N  bodies  within  an  area  which  is  called  the  searching  space  or  working  space.  In  the 
definition  of  the  searching  space,  we  use  the  following  a  priori  assumption:  during  the  time 
interval  of  interest,  the  bodies  in  question  never  leave  the  searching  space,  so  its  area  remains 
unchanged  throughout  the  evolution  of  the  problem. 

Let  us  consider  the  set  of  surface  nodes  corresponding  to  the  surfaces  of  the  N  bodies  and 
denote  by  Utotai  the  size  of  this  set  of  nodes. 

Remark  7.1  Prom  the  point  of  view  of  the  previous  chapters,  this  set  of  nodes  corresponds 
to  the  surface  nodes  of  the  slave  and  master  bodies  which  participate  in  the  contact  problem 
prior  to  any  closest  point  projection  procedure.  The  set  of  surface  nodes  of  the  master  body 
must  not  be  confused  with  the  master  segment  identified  by  the  closest  point  projection. 

NOTE:  In  the  rest  of  the  chapter,  we  give  the  name  particles  to  the  set  of  Utotai  nodes 
belonging  to  the  surfaces  of  the  bodies  participating  in  the  contact  problem,  and  we  give  the 
name  nodes  to  the  nodes  of  the  binary  tree  using  in  storing  our  data. 

The  root  of  the  binary  tree  is  associated  with  our  chosen  initial  searching  space.  Let  us 
consider  a  particle  with  its  corresponding  current  coordinates.  We  divide  the  working  space 
in  half  in  a  particular  direction,  say  vertically,  and  then  determine  on  which  side  this  particle 
falls;  then,  we  create  a  child  which  we  associate  with  that  half  space.  We  think  of  this 
procedure  as  inserting  the  particle  in  a  node.  Next,  we  consider  another  particle.  If  this  new 
particle  falls  on  the  same  half  space  as  the  previous  one,  we  subdivide  the  half  space  and 
we  move  the  previously  inserted  particle  into  the  corresponding  quarter  space,  and  insert 
the  new  particle  into  the  quarter  space  where  it  belongs.  The  procedure  is  repeated  until 
each  particle  resides  on  a  leaf.  On  the  other  hand,  if  the  new  particle  falls  on  the  opposite 
half  space  as  the  previous  one,  we  just  create  the  second  child  in  the  tree  and  continue  the 
process  by  considering  the  next  particle  in  the  queue. 

Remark  7.2  Notice  that  we  do  not  create  a  child  node  in  our  tree  unless  there  is  a  particle 
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Figure  7.2:  Schematic  drawing  of  a  set  of  particles  within  a  square  working  space 

occupying  it.  This  is  a  modified  way  of  building  a  tree  and  avoids  the  creation  of  nodes  which 
later  would  have  to  be  deleted  because  they  are  empty.  Other  authors  use  different  versions 
of  tree  building  algorithms  which  are  more  suitable  for  their  purposes  (see  e.g.  Bonet  & 
Peraire[6]). 

When  we  finish  inserting  all  the  particles  in  the  tree,  they  all  reside  on  leaves,  i.e.  there  is  a 
unique  particle  residing  in  each  cell. 

7.3.2. 1  Construction  of  the  binary  tree  example 

In  Figures  7.2  and  7.3  we  show  an  example  of  a  set  of  particles  which  have  been  inserted 
into  a  tree. 
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Figure  7.3:  Binary  tree  structure  corresponding  to  the  example  shown  in  Figure  7.2 

7.3. 2. 2  Binary  tree  construction  program 

Given  a  list  of  particles  coordinates,  the  binary  tree  can  be  constructed  recursively  as  follows: 

Procedure  BintreeBuild 
Bintree  =  {empty} 
do  i  1 .  ^total 

Call  BinlnsertCi  ,root) 
end  do 

Procedure  Binlnsert(i,n)... insert  particle  i  in  node  n 
if  the  subtree  rooted  at  n  contains  more  than  1  'particle 
Determine  which  child  c  of  node  n  particle  i  lies  in 
Call  Binlnsert(i,c) 

else  if  the  subtree  rooted  at  n  contains  1  particle,  n  is  a  leaf 
Consider  n^s  2  children,  create  the  child  in  which 
the  particle  already  in  n  liesand  move  particule  i  into  it. 

Let  c  be  the  child  in  which  peirticle  i  lies 
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Call  Binlnsert(i,c) 

else  if  the  subtree  rooted  at  n  is  empty,  n  is  a  leaf 
Store  particle  i  in  node  n 
end  if 

7.4  Sorting  using  binary  tree 

In  this  section,  we  explain  how  best  to  use  the  data  stored  in  the  binary  tree  in  a  general 
sorting  algorithm  (see  Munjiza  et  Al[24]).  Once  the  binary  tree  has  been  constructed,  as 
described  in  the  previous  sections,  we  proceed  to  the  actual  sorting  phase  of  the  algorithm. 
We  search  for  possible  contact  pairs  between  any  of  the  Utotai  particles  and  the  N  bodies. 
Recall  that  the  particles  are  the  actual  surface  nodes  (in  the  finite  element  sense)  of  these 
bodies. 

Remark  7.3  We  assume  that  a  node  cannot  penetrate  the  body  it  belongs  to,  thus  elimi¬ 
nating  self  contact  cases  from  the  algorithm. 

To  check  for  contact  at  this  stage,  we  construct  a  buffer  zone  around  each  body  to  simplify 
the  search  for  contact  [24].  An  example  of  a  buffer  zone  is  the  circumscribing  rectangle 
which  is  aligned  with  the  reference  coordinate  system  (see  Figure  7.4).  The  coordinates  of 
this  rectangle  are  the  minimum  and  maximum  coordinates  in  each  space  dimension  that 
bound  the  object.  In  our  work,  we  have  used  this  type  of  buffer  zone,  which  does  not  assume 
any  a  priori  knowledge  of  the  problem. 

Remark  7.4  The  calculation  of  the  buffer  zone  is  related  to  the  sorting  scheme,  since  it 
depends  on  how  often  one  performs  a  sorting  procedure.  If  a  sorting  procedure  is  performed 
every  n  time  steps,  the  buffer  zone  is  chosen  so  that  its  thickness  b  is  at  least  b  >  nvmax^t, 
where  Vmax  is  the  maximum  velocity  among  all  the  bodies. 
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Figure  7.4:  Schematic  drawing  of  a  body  and  its  corresponding  buffer  zone 

7.4.1  Traversing  the  binary  tree 

The  contact  search  is  performed  by  traversing  the  tree.  The  binary  tree  is  traversed  N  times 
(once  for  each  body).  We  visit  each  node  on  the  tree,  and  it  has  associated  with  it  a  square  or 
rectangular  cell.  The  target  body  with  its  designated  buffer  zone  is  checked  for  superposition 
with  this  cell.  Efficiency  is  improved  if  we  take  into  account  that  if  no  superposition  with  a 
node  is  detected,  then  no  further  contact  detection  need  be  performed  for  the  corresponding 
subtree  rooted  at  that  node. 

When  we  reach  a  leaf  of  the  tree  which  contains  an  actual  particle  (slave  or  master  node), 
the  particle  is  checked  for  inclusion  within  the  body  in  question. 

To  traverse  the  tree  we  use  the  recursive  preorder  scheme.  (See  [17]  for  the  basics  on 
traversing  trees  and  preorder  schemes.) 

7.4.1. 1  Preorder  scheme 

The  pseudo  code  is  presented  below 

subroutine  preorder  (current-address) 
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call  visit  (Visit  node  stored  at  current-address) 
if  (cent  =  .true.)  then 
if  (current-left-child  exists)  then 
call  preorder (current-left-child-address) 
end  if 

if  (current-right-child  exists)  then 
call  preorder (current-right-child-address) 
end  if 
end  if 
stop 
end 

Visiting  a  node  means  checking  for  superposition  of  two  rectangles  (the  subpartition  as¬ 
sociated  with  the  tree  node  and  the  buffer  zone  rectangle)  or  checking  whether  a  particle 
lies  within  a  rectangle  (the  finite  element  node  and  the  buffer  zone  rectangle).  During  the 
‘Visit” ,  then,  we  recognize  a  pair  consisting  of  a  particle  and  a  body  as  being  in  possible 
contact.  This  pair,  which  we  call  bp-pair  for  body-particle-pair,  is  then  added  to  the  list  of 
possible  pairs  in  contact. 

Remcirk  7.5  Even  though  we  do  not  assign  slave  and  master  categories  to  the  particles 
in  our  problem  before  the  sorting  phase,  the  categorization  comes  up  automatically.  The 
closest  point  projection  phase  of  the  algorithm  will  consider  that  the  particle  in  the  bp-pair 
is  the  slave  node  to  be  checked  for  penetration  into  the  body  belonging  to  the  bp-pair. 


7.5  Implementation  of  a  sorting  algorithm 

The  type  of  database  involved  in  the  sort/search  procedure  requires  the  use  of  data  structures 
and  pointers.  FORTRAN  90  provides  us  with  these  necessary  tools. 

The  design  and  definition  of  our  database  structure  are  of  importance,  as  they  affect  the 
eflaciency  and  memory  use  of  the  algorithm.  Below,  we  present  some  of  the  most  important 
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structures  and  issues  involved  in  the  implementation. 

7.5.1  Binary  tree 

A  node  in  the  tree  is  defined  as  follows: 
type  node 

type  (particle)  .pointer  ::  part-in-node 
real  (8),  dimens  ion  (2)  ::  xmin,xmax 
type  (node)  .pointer  : :  parent, left-child, right-child 
end  type  node 

Pointers  are  used  to  point  at  other  pieces  of  data,  thus  eliminating  the  need  to  make  multiple 
copies  of  the  data.  For  example,  when  traversing  the  binary  tree,  we  need  to  go  down  to 
the  children  of  a  certain  node,  so  we  include  in  the  node  data  structure  two  pointers,  one  to 
each  child. 

The  first  part  of  the  node  data  structure  is  the  part-in-node,  which  is  a  pointer  to  a  particle 
type  structure.  If  the  node  is  not  a  leaf  then  this  pointer  is  nullified  (does  not  point  to  any 
data).  When  this  happens,  the  pointer  is  said  to  be  disassociated. 

7.5.2  Body,  surface  and  particle 

A  body  is  formed  by  a  set  of  surfaces  and  each  one  of  these  surfaces  consists  of  a  set  of 
particles.  Following  this  reasoning,  we  can  think  of  a  particle  as  being  a  substructure  of  a 
surface  and  a  surface  as  a  substructure  of  a  body. 

Keeping  this  in  mind,  we  define  the  following  structures: 
type  body 

integer  : :  num-surf 

type  (surface)  . dimension  ( : ) . pointer  : :  id- surf 
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real (8),  dimension(2)  ::  lim-min,lim-max 
end  type  body 

The  integer  num-surf  indicates  the  number  of  surfaces  that  define  the  body.  The  real 
numbers  lim-min  and  lim-max  contain  the  data  of  the  buffer  zone  rectangle.  The  pointer 
array  id-surf  points  to  each  surface  data  structure,  as  described. 

Remark  7.6  The  capability  of  defining  a  body  by  means  of  multiple  surfaces,  enables  the 
scheme  to  simulate  a  body  with  different  surface  properties.  For  example,  a  body  may  have 
different  surface  finishes,  each  with  a  different  friction  coefficient. 

type  surface 

integer  : :  num-part 

type  (surface)  ,dimension(:)  .pointer  ::  id-part 
end  type  surface 

The  integer  num-part  indicates  the  number  of  particles  that  define  the  surface.  The  pointer 
array  id-part  points  to  each  particle  data  structure  as  described. 

type  particle 

integer  : :  glob-num,b-id,s-id 
real (8)  x,y 
end  type  particle 

The  integers  glob-num,  b-id  and  s-id  represent  the  global  node  number  (in  the  finite  element 
sense),  the  body  number  and  surface  number  which  this  particle  belongs  to,  respectively. 

7.5.3  Linked  lists 

The  output  of  the  sorting  phase,  the  set  of  bp-pairs,  is  stored  in  a  linked  list.  The  data 
structure  bp-pair  is  defined  as  follows: 
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type  bp-pair 

type  (body)  .pointer  ::  bp-pair-body 
type  (particle)  .pointer  bp-pair-part 
type  (bp-history)  .pointer  ::  bp-pair-histl,  bp-pair-hist2 
end  type  bp-pair 

The  pointer  bp-pair-body  points  to  the  body  data  that  is  in  contact  with  the  particle  whose 
data  is  pointed  to  in  turn  by  bp-pair-part.  The  data  bp-pair-histl  and  bp-pair-hist2  cor¬ 
respond  to  two  history  arrays  belonging  to  data  at  tn  and  tn+i,  respectively,  and  have  the 
following  data  structure: 

type  bp-history 
logical  : :  cont 
integer. dimension (3)  : :  ixl 
real (8)  ::  norin(2),  gapd 

. . .  other  data  needed  for  the  calculation  of  the  force  and  tangent  matrix  contribution 
end  t3rpe  bp-history 

The  logical  variable  cont  is  the  contact  flag  for  the  considered  time  step  considered.  The 
integer  array  ixl  contains  the  global  node  numbers  (in  the  finite  element  sense)  of  the  three 
nodes  that  form  our  three-node  contact  element.  The  real  variable  gapd  is  the  dynamic  gap 
which  we  need  to  keep  track  of  for  our  restoring  scheme. 

When  the  sorting  scheme  detects  a  possible  contact,  the  bp-pair  is  added  to  a  linked  list. 
When  no  contact  is  detected  by  the  closest  point  projection  procedure,  these  bp-pairs  are 
taken  out  from  the  linked  list  and  the  space  allocated  for  them  is  then  deallocated.  This 
procedure  has  the  advantage  of  using  memory  in  a  dynamic  way  thus  saving  memory  storage. 

A  linked  list  is  composed  of  list-nodes.  The  structure  of  a  list-node  is  as  follows: 
type  list-node 

type  (bp-pair)  .pointer  ::  elem 
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type  (list- node) , pointer  ::  prev,next 
end  type  list-node 

The  pointer  elem  points  to  the  information  of  the  bp-pair  that  is  contained  in  the  given 
list-node  and  the  pointers  prev  and  next  point  to  the  previous  and  next  list-nodes  in  the 
linked  list,  respectively. 

7.5.4  Auxiliary  tools 

In  addition  to  all  the  previously  defined  data  structures,  we  use  some  auxiliary  arrays.  One 
of  the  most  important  is  an  array  of  pointers  called  pairs-table,  of  dimension  ntotai  x  N. 
These  pointers  are  initially  nullified,  but  if  a  particular  bp-pair  (with  particle  number  np 
and  body  number  nb)  ,  was  active  at  the  converged  state  of  the  previous  time  step,  then 
pairs-table(n6,np)  points  to  that  particular  bp-pair,  and  the  algorithm  neither  allocates 
more  memory  nor  generates  a  second  copy  of  the  same  bp-pair.  On  the  other  hand,  if  pairs- 
table(n&,np)  is  nullified  then  it  is  evident  that  this  particular  bp-pair  was  not  active  at  the 
previous  time  step,  and  the  algorithm  must  allocate  memory  to  create  it. 

7.5.5  Basic  algorithm 

With  the  salient  points  of  our  database  explained,  we  now  detail  the  basic  sorting/search 
algorithm  used  in  the  implementation  of  our  contact/impact  scheme. 

In  the  box  above,  we  denote  by  a%b  the  element  6  of  a  certain  type,  where  a  is  of  that  type. 

Example:  Using  FORTRAN  90,  we  may  declare 

type{node)  ::  a  (o  is  a  node  in  the  binary  tree)  , 
so  that  we  have  access  to  the  left  or  right  child  of  a,  respectively  by 

I  =  a%left  —  child  or  r  =  a%right  —  child 
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BASIC  ALGORITHM  FOR  THE  SORT/SEARCH  PROCEDURE 

•  Create  root 

•  Update  particle  coordinate 

•  Calculate  buffer  boxes 
DO  z  =  1  :  ntotai 

•  Insert  particle  i  in  the  binary  tree 
END  DO 

BO  i  =  l  :N 

•  Traverse  the  tree  and  check  for  contact  between  body  i  and 
the  particles  in  each  leaf  of  the  tree 

END  DO 

•  Delete  binary  tree 

DO  WHILE  list%bp-pair  exists 

IF  (list%bp-pair  existed  at  time  (check  pairs-table))  THEN 

•  Perform  closest  point  projection  at  and 
calculate  dynamic  gap  Check  for  contact 

ELSE 

•  Perform  closest  point  projection  at  and 
calculate  real  gap  Qn 

•  Perform  closest  point  projection  at  and 
calculate  dynamic  gap  (with  p„).  Check  for  contact 

END  IF 

IF  (cont  =  .true.)  THEN 

•  Calculate  the  force 
END  IF 

list  list%next 
END  DO 


Table  7.1:  Basic  algorithm  for  the  sort /search  procedure 
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where  I  and  r  are  declared  as 

type{node)  ::  l,r  [I  and  r  are  nodes  in  the  binary  tree)  , 

Remark  7.7  Notice  that  the  binary  tree  is  deleted  as  soon  as  we  form  the  linked  list  con¬ 
taining  all  the  bp-pairs  in  possible  contact. 


7.6  Computation  times 

As  stated  above,  the  addition  of  a  sorting  phase  decreases  the  computational  effort  involved 
in  the  contact  detection  part  of  the  contact  algorithm  in  multibody  problems. 

In  this  section,  we  measure  (1)  the  average  CPU  time  for  sorting/searching  in  each  time 
step,  and  (2)  the  average  CPU  time  for  the  entire  step,  for  a  sample  of  problem  while  we 
vary  the  number  of  bodies  involved. 

These  measurements  give  us  an  estimate  of  the  speed  increase  provided  by  the  sorting  pro¬ 
cedure  as  the  number  of  bodies  increases.  It  is  also  of  interest  to  see  the  fraction  of  time  the 
contact  detection  part  takes  within  an  overall  time  step. 

Theoretical  calculations  show  that  for  a  homogeneous  concentration  of  particles,  which 
should  produce  a  well-balanced  tree,  the  times  for  the  building  and  traversing  the  tree 
would  behave  as  0{N\og{N)),  where  N  is  the  number  of  particles.  The  detection  part  of 
the  algorithm  with  no  sorting  phase,  i.e.  an  exhaustive  search,  behaves  as  0{N^).  In  Table 
7.2,  we  show  the  CPU  time  in  seconds  for  different  problems. 

A  regression  analysis  has  been  performed  with  the  times  from  Table  7.2  to  evaluate  the 
behaviour  of  each  algorithm.  With  the  sorting  algorithm,  the  tested  times  for  detetction 
followed  the  expected  logarithmic  dependence  on  N  with  a  correlation  coefficient  oi  R  — 
0.99985,  while  the  times  reuired  wihtout  sorting  matched  a  quadratic  curve  with  a  correlation 
coefficient  ol  R  —  0.999991.  Figures  7.5  and  7.6  show  the  curves  obtained  through  regression 
analysis  of  the  data  in  Table  7.2. 

The  ratio  between  the  CPU  time  for  the  detection  phase  and  the  CPU  time  for  the  solver 
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N 

Sorting 

No  sorting 

detection 

CPUio,,, 

detection 

CPUtota/ 

13 

0.0085 

0.07 

0.1244 

0.29 

53 

0.0361 

0.90 

2.7050 

3.60 

104 

0.0732 

1.80 

10.8600 

12.80 

196 

0.1386 

3.60 

40.3800 

44.23 

400 

0.3030 

7.64 

900 

0.7540 

16.60 

- 

- 

Table  7.2:  Computing  times  in  seconds  for  various  numbers  of  bodies. 


Figure  7.5:  Contact  detection  CPU  time  for  the  algorithm  with  sorting  procedure.  Compu¬ 
tational  data  (  o  );  regression  analysis  ( — ). 


7.6.  COMPUTATION  TIMES 


80  100  120  140  160  180  200 

number  of  bodies 


Figure  7.6:  Contact  detection  CPU  time  for  the  algorithm  with  no  sorting  procedure.  Com¬ 
putational  data  (  o  );  regression  analysis  ( — ). 
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N 

•^sorting 

^no  sorting 

13 

0.0012 

0.74 

53 

0.0418 

2.91 

104 

0.0424 

5.60 

196 

0.0400 

10.49 

400 

0.0413 

- 

900 

0.0476 

- 

Table  7.3:  Computing  times  in  seconds  for  various  numbers  of  bodies. 

procedure  is  also  of  interest,  since  without  the  addition  of  a  sorting  phase  the  detection 
algorithm  would  dominate  the  CPU  time  for  large  N.  We  denote  this  ratio  by  A,  i.e., 


A 


case 


CPU 


^  detection _ 

case  _ ppTTcase 

total  ^detection 


(7.1) 


where  case  refers  to  sorting  or  no  sorting.  In  Table  7.3  we  show  these  values  for  different 
number  of  bodies. 


One  may  observe  that  the  ratio  A  increases  linearly  with  N  for  the  non-sorting  algorithm, 
while  the  value  is  independent  of  N  when  sorting  is  added  to  the  search  algorithm. 
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Chapter  8 


Numerical  Examples 


8.1  Introduction 

In  this  section  we  present  a  variety  of  numerical  simulations  with  which  we  evaluate  the 
performance  of  the  various  schemes  presented  in  this  work.  The  three  sections  detail  a  variety 
of  benchmark  problems  to  test  the  correctness  of  the  solutions  and  some  other  problems  to 
compare  the  quality  of  the  solutions  with  other  more  conventional  contact/impact  schemes. 
The  last  section  primarily  focuses  on  the  versatility  of  the  sort /search  scheme  of  the  contact 
detection  part  of  the  algorithm  in  cases  of  multibody  problems.  The  algorithms  formulated 
in  this  work  have  been  implemented  in  the  general  finite  element  code  FEAP  developed  by 
Professor  Robert  L.  Taylor  of  the  University  of  California  at  Berkeley. 


8.2  Prict ionless  contact  simulations 

8.2.1  Impact  of  a  rod  on  a  rigid  wall 

The  purpose  of  this  simulation  is  to  show  the  important  role  that  the  energy  restoring 
property  of  our  scheme  has  on  the  quality  of  the  solution  and  on  the  overall  stability  of  the 
simulation.  The  contact  problem  of  elastic  bodies  may  be  thought  of  as  consisting  of  two 
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Figure  8.1:  Impact  of  an  elastic  rod  against  a  rigid  wall.  Schematic  drawing  of  the  problem. 

parts:  the  continuum  part,  which  deals  with  the  elastic  deformation,  and  the  contact  part, 
which  deals  with  the  unilateral  constraint  (impenetrability). 

We  solve  the  continuum  problem  using  the  conservative  scheme  presented  in  Chapter  3  and 
developed  by  SiMO  &  Tarnow  [30],  and  compare  these  results  with  the  ones  obtained 
through  HHT  scheme,  also  known  as  the  a-method,  presented  in  Hilber  et  al  [13].  When 
solving  the  contact  part  of  the  problem  we  use  the  restoring  scheme  presented  in  Chapter  5 
and  the  mid-point  contact  scheme  developed  by  Laursen  &  SiMO  [20]. 

The  HHT  scheme  is  dissipative  for  linear  problems,  hence  unconditionally  stable,  but  it  this 
dissipative  property  for  general  nonlinear  problems.  Thus,  in  the  case  of  a  contact  problem, 
one  may  find  an  increase  in  energy  if  the  scheme  is  not  used  in  combination  with  an  adequate 
contact  algorithm  to  deal  with  the  nonlinearity  inherent  to  the  contact  constraint. 

Let  us  consider  the  one  dimensional  model  of  a  rod  impacting  a  rigid  wall  which  can  be 
solved  using  various  combinations  of  continuum  and  contact  algorithms. 

The  problem  is  sketched  in  Figure  8.1.  The  rod  used  in  the  simulations  has  length  L  =  1.0 
and  is  discretized  with  100  isoparametric  elements.  The  Young’s  modulus  is  £■  =  1.0  and 
the  density  is  p  =  1.0.  The  initial  velocity  before  impact  is  vq  =  —0.5.  The  rod  is  initially 
at  a  distance  of  7.5  x  10~^  from  the  wall.  The  wave  velocity  within  the  bar  is  c  = 
and  the  simulation  is  performed  using  a  discrete  CFL  =  ^  =  2,  where  h  is  the  element 
length  in  the  finite  element  discretization.  Analytical  calculations  show  that  the  rod  will 
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remain  in  contact  with  the  rigid  wall  for  a  period  of  time  T  =  ^,  so  that  in  our  problem 
we  have  T  =  2.  The  simulations  should  then  yield  g  =  0  and  g  =  0  for  t  G  [e,  e  +  2], 
where  e  =  — —  (the  time  the  rod  takes  to  reach  the  wall  from  its  initial  position).  The 
simulations  are  performed  with  a  time  step  At  =  0.0025;  since  T  >  At,  the  imposition  of 
the  gap  velocity  constraint  becomes  necessary. 

We  assume  a  linear  model  so  that  we  only  introduce  the  nonlinearity  arising  from  the  contact 
phenomenon.  We  consider  the  three  parameter  family  of  time  discretizations  (see  Hilber 
ET  AL  [13])  to  accomodate  dissipative  schemes  when  we  deal  with  the  continuum  problem. 
Thus,  we  have 


/c,n+|  =  Man+i  +  IC[adn+i  +  (1  ~  Ci)dn]  ,  (8.1) 

dn+1  =  dn  +  Atvn  +  - A<^[2/?a„+i  +  (1  —  2/5)a„]  ,  (8.2) 

Vn+1  =  Vn  +  At[7a„+1  +  (1  -  7)a„]  ,  (8.3) 

where  M  denotes  the  standard  mass  matrix  and  K  denotes  the  usual  stiffness  matrix  of 
linear  elasticity.  The  equilibrium  equation  8.1  is  written  in  the  form  presented  in  SiMO  ET 
AL  [31],  where  the  parameter  a  in  [13]  is  equivalent  to  a  —  1  in  [31]. 

We  also  consider  a  standard  mid-point  contact  scheme  by  imposing  the  contact  constraint 
at  tn+a,  where  a  is  consistent  with  to  the  one  used  in  equation  8.1.  The  contact  pressure  t^ 
then  takes  the  following  form: 


tjV  —  I^N9n+a  )  (^-4) 

where  gn+a  is  the  (real)  gap  calculated  using  a  closest  point  projection  in  the  configuration 
at  time  tn+a- 

We  consider  the  following  combination  of  schemes: 

Algorithm  1  Trapezoidal  rule:  a  =  1.0,  /?  =  0.25  and  7  =  0.5  for  the  continuum  part;  and 
the  standard  penalty  scheme  described  above  for  the  contact  part. 
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Algorithm  2  Mid-point  rule:  a  —  0.5,  ^  =  0.5  and  7  =  1.0  for  the  continuum  part;  and 
the  standard  penalty  scheme  described  above  for  the  contact  part. 

Algorithm  3  HHT:  a  =  0.51,  /?  =  0.555025  and  7  =  0.99  for  the  continuum  part;  and  the 
standard  penalty  scheme  described  above  for  the  contact  part. 

Algorithm  4  Mid-point  rule:  a  =  0.5,  =  0.5  and  7  =  1.0  for  the  continuum  part;  and 

the  new  energy  restoring  contact  scheme  for  the  contact  part. 

Algorithm  5  Mid-point  rule:  a  =  0.5,  13  =  0.5  and  7  =  1.0  for  the  continuum  part;  and 
the  new  dissipative  contact  scheme  for  the  contact  part. 

Note  that  the  conserving  elastic  algorithm  described  in  Chapter  3  yields  the  mid-point  and 
trapezoidal  rules  when  applied  to  a  linear  elastic  problem,  so  these  are  energy  conserving 
schemes.  On  the  other  hand,  the  HHT  scheme  is  energy  dissipative  for  linear  problems. 
As  we  illustrate  with  this  numerical  simulation.  Algorithms  1,  2  and  3  lose  the  energy 
conserving  property  when  the  contact  nonlinearity  is  added  to  the  problem  and  treated 
using  the  standard  mid-point  contact  penalty  regularization.  The  results  also  show  that 
when  either  of  the  continuum  schemes  is  combined  with  the  appropriate  contact  scheme,  the 
energy  conserving  property  remains  valid. 

Figure  8.2.1  shows  the  results  for  the  trapezoidal  rule  (Algorithm  1)  compared  with  the 
proposed  scheme  (Algorithm  4).  While  Algorithm  4  successfully  constrains  the  gap 
and  the  gap  velocity  during  the  time  of  contact,  the  trapezoidal  rule  evidences  a  marked 
oscillatory  behaviour.  This  type  of  behaviour  is  also  reflected  in  the  pressure  plot  while 
Algorithm  4  yields  a  good  approximation  of  the  step  function.  The  energy  predicted  by 
Algorithm  1  grows  almost  monotonically  during  the  time  of  contact,  while  the  proposed 
scheme  (Algorithm  4)  restores  the  energy  after  the  release  of  the  rod. 

As  one  may  observe  in  Figure  8.2.1,  the  oscillatory  behaviour  is  repeated  by  the  standard 
contact  scheme  (Algorithm  2).  The  energy  grows  when  the  mid  point  rule  (energy  conserv¬ 
ing)  is  using  in  conjunction  with  the  standard  penalty  scheme  instead  of  the  energy  restoring 
one. 


93 


8.2.  FRICTIONLESS  CONTACT  SIMULATIONS 


Time 


Time 


Figure  8.5:  Impact  of  an  elastic  rod  against  a  rigid  wall.  Plots  comparing 
Algorithm  5  with  ^  =  0.9  (•  •  • )  and  Algorithm  4  ( — ). 

For  the  purpose  of  keeping  the  oscillations  under  control  when  using  the  standard  penalty 
scheme,  one  can  use  a  dissipative  scheme  to  deal  with  the  continuum  part  as  in  the  case  of 
Algorithm  3  (see  Figure  8.4).  The  aim  of  this  comparison  is  to  prove  that  the  dissipation 
provided  by  the  HHT  scheme  is  not  enough  to  prevent  the  growth  of  the  total  energy  during 
the  time  of  contact.  As  observed  in  Figure  8.4,  the  combination  of  schemes  successfully 
eliminates  the  oscillations  in  the  evolution  of  the  gap  function,  but  they  are  still  present  in 
the  gap  velocity  and  pressure  plots.  The  energy  increases  during  a  short  interval  of  time, 
confirming  our  hypothesis. 

The  last  example,  shown  in  Figure  8.5,  demonstrates  the  effect  of  using  a  dissipative  scheme 
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to  deal  with  an  impact  case  which  has  high  frequency  energy  modes.  The  dissipation  ef¬ 
fectively  eliminates  even  the  small  oscillations  present  in  the  energy  restoring  scheme.  The 
possibility  of  introducing  dissipation  in  a  regulated  way  by  changing  the  value  of  the  param¬ 
eter  0  makes  the  algorithm  a  very  attractive  choice  for  problems  where  the  high  frequency 
energy  modes  tend  to  overwhelm  the  results  of  the  simulation. 


8.2.2  Impact  of  two  nonlinear  elastic  cylinders 

Consider  the  impact  if  two  nonlinear  elastic  cylinders  in  plane  strain.  The  cylinders  have  a 
diameter  of  2.0  and  are  discretized  using  bilinear  quadrilateral  finite  elements.  The  Saint- 
Venant  Kirchhoff  material  model  is  assumed  for  both  cylinders.  We  simulate  the  impact  of 
these  two  cylinders  with  two  very  different  sets  of  material  properties.  In  this  way,  the  ability 
of  the  scheme  developed  in  this  thesis  to  simulate  very  different  time  scales  and  degrees  of 
deformation  is  illustrated. 

8.2.2.1  Quasi-rigid  cylinders 

Both  cylinders  have  Lame  constants  A  =  2.0  x  10^  and  =  1.0  x  10^,  and  density  p  =  1.0. 
These  properties  make  the  cylinders  behave  in  a  quasi-rigid  fashion  so  that  the  duration  of 
contact  is  very  small  compared  to  the  overall  period  of  simulation.  Following  the  criteria 
described  in  Chapter  5,  the  satisfaction  of  the  gap  velocity  constraint  becomes  quite  unnec¬ 
essary.  Thus,  the  penalty  parameters  used  are  as  follows:  kjv  =  1-0  x  10®  and  rrip  =  0.0  (no 
satisfaction  of  the  gap  velocity  constraint).  A  constant  time  step  At  =  0.1  is  used  in  the 
simulation.  The  left  cylinder  is  given  an  initial  velocity  vq  =  (1.0, 0.2);  where  thereafter  a 
series  of  collisions  among  the  cylinders  and  the  surrounding  walls  may  be  observed  in  Figure 
8.6. 

The  aim  of  this  example  is  to  show  that  the  unchecked  growth  of  energy  affects  the  dynamics 
of  the  two  cylinders  so  that  after  the  first  impact  the  evolution  of  the  system  depends 
significantly  on  which  algorithm  is  used. 
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Figure  8.6;  Impact  of  quasi-rigid  cylinders.  Evolution  of  the  system 
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Figure  8.7:  Impact  of  quasi-rigid  cylinders.  Evolution  of  the  linear  momenta,  angular  mo¬ 
mentum  and  the  energy  of  the  system.  Plots  omparing  Algorithm  2  (•  •  • )  and  the  proposed 
scheme  ( — ). 
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8.7  W6  plot  tho  linosr  inornonta.,  angular  momontum  and  tli6  onorgy  6 volution  for 
the  system,  using  the  proposed  scheme  and  Algorithm  2  (which  uses  the  mid-point  rule 
to  impose  the  contact  constraint).  One  may  observe  that  the  first  impact  occurs  between 
the  two  cylinders  at  approximately  t  =  1.8,  when  the  results  show  a  jump  in  the  energy 
evolution  for  the  mid-point  rule.  The  proposed  scheme  shows  a  small  decrease  in  energy 
during  the  contact  period  before  the  energy  is  restored  to  its  original  value.  The  second 
impact  of  the  problem  occurs  at  approximately  t  =  2.8  between  one  cylinder  and  a  rigid 
wall.  Theoretical  results  predict  that  the  linear  momentum  in  the  x  direction  should  change 
sign  with  a  small  variation  in  its  magnitude  while  the  linear  momentum  in  the  y  direction 
should  remain  unchanged.  In  the  linear  momentum  plots  we  observe  that  this  behaviour  is 
captured  by  our  scheme.  With  the  mid-point  rule,  the  linear  momentum  in  the  x  direction 
changes  sign  but  its  value  is  almost  doubled.  At  this  instant,  the  energy  shows  a  marked 
increase. 

The  trajectory  of  the  cylinders  with  the  standard  mid-point  rule  scheme  is  distinctly  different 
from  the  one  shown  in  Figure  8.6,  so  the  instants  at  which  impact  occurs  are  also  quite 
different,  as  evidenced  by  the  plots  in  Figure  8.7. 

8. 2.2.2  Soft  cylinders 

The  second  set  of  material  parameters  for  both  cylinders  are  as  follows:  Lame  constants 
A  =  130.0  and  /i  =  43.33,  and  density  p  =  8.93.  With  this  set  of  values,  the  cylinders  will 
undergo  large  deformations  and  a  long  contact  duration.  The  left  cylinder  is  given  an  initial 
velocity  =  (1.0, 0.1). 

The  evolution  of  the  system  is  shown  in  Figure  8.8.  One  may  observe  that  the  time  interval  of 
interest  in  this  case  is  of  the  same  order  as  the  contact  time.  Also,  the  amount  of  deformation 
;of  the  cylinders  is  significant. 

Figure  8.9  shows  a  significant  increase  of  energy  during  the  contact  time  when  predicted  by 
the  mid-point  scheme,  whereas  the  proposed  scheme  restores  the  energy  to  its  initial  value 
at  the  end  of  the  contact  interval. 
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Figure  8.8:  Impact  between  two  soft  elastic  cylinders.  Evolution  of  the  dynamical  system. 
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8.3  Frictional  contact  simulations 

8.3.1  Impact  of  a  cylinder  against  a  rigid  wall 

We  present  in  this  section  the  results  obtained  by  means  of  the  proposed  scheme  for  the 
problem  of  an  elastic  cylinder  impacting  a  rigid  wall  at  an  angle.  The  problem  has  been 
solved  using  both  frictionless  and  frictional  conditions. 

We  have  used  a  Saint- Venant  Kirchhoff  model  with  Lame  constants  A  =  130.0  and  G  =  43.33, 
and  a  density  p  =  8.93.  The  cylinder  has  a  radius  R  =  1.0  and  is  given  an  initial  velocity 
VO  =  (0.4,  -0.4). 

In  the  frictional  case  the  friction  constant  p,  =  0.2  has  been  used,  with  the  following  penalty 
parameters:  =  kt  =  10'^  with  0  =  0.5.  The  continuum  contributions  of  the  problem  were 

solved  using  the  energy  conserving  scheme  detailed  in  [30]. 

To  compare  the  overall  performance  of  our  scheme,  we  have  also  solved  the  problem  using  the 
standard  mid-point  rule  for  the  contact  and  friction  contributions  and  the  energy  conserving 
scheme  for  the  continuum  contributions. 

Figures  8.10  and  8.11  show  the  results  of  the  proposed  scheme  for  three  instances  of  the  fric¬ 
tionless  and  frictional  problem,  respectively.  Notice  the  absence  of  rotation  in  the  frictionless 
case. 

We  observe  that  the  mid-point  rule  yields  an  increase  of  energy,  as  shown  in  Figure  8.12  for 
the  frictionless  case.  An  increase  is  predicted  even  when  there  is  friction  involved  (Figure 
8.13),  which  should  add  dissipation  to  the  system.  Therefore,  we  conclude  that  the  use  of 
the  dynamic  contact  quantities,  for  both  the  frictionless  and  the  frictional  case,  contribute 
to  the  overall  stability  of  the  contact  scheme. 
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Figure  8.10:  Impact  of  an  elastic  cylinder  against  a  rigid  wall.  Three  instances  of  the 
evolution  of  the  frictionless  case,  at  times  t  =  0, 6, 12. 


Figure  8.11:  Impact  of  an  elastic  cylinder  against  a  rigid  wall.  Three  instances  of  the 
evolution  of  the  frictional  case  {fj.  =  0.2),  at  times  t  =  0, 6, 12. 
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Figure  8.12:  Impact  of  an  elastic  cylinder  against  a  rigid  wall.  Energy  evolution  of  the 
frictionless  case.  Plots  comparing  the  proposed  scheme  ( — )  and  the  mid-point  rule  (•••)• 


Figure  8.13:  Impact  of  an  elastic  cylinder  against  a  rigid  wall.  Energy  evolution  for  the 
frictional  case  (//  =  0.2).  Plots  comparing  the  proposed  scheme  ( — )  and  the  mid-point  rule 

(•••)• 
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Figure  8.14:  Forging  of  an  elastic  block  against  a  rigid  foundation.  Schematic  drawing  of 
the  problem. 

8.3.2  Forging  of  an  elastic  block  against  a  rigid  foundation 

The  purpose  of  this  simulation  is  to  show  the  performance  of  our  scheme  for  a  frictional 
problem  and  to  compare  the  results  with  other  frictional  schemes.  In  addition,  this  example 
shows  that  the  proposed  scheme  can  be  used  not  only  in  dynamic  problems  but  also  in 
quasi-static  problems. 

Consider  the  problem  of  an  elastic  block  pressed  against  a  rigid  foundation.  The  block  is 
pulled  by  a  tangential  force  uniformly  distributed  along  one  of  the  sides  of  the  block,  as 
shown  in  Figure  8.14. 

The  problem  has  been  solved  using  bilinear  quadrilateral  elements  within  the  linear  elastic 
formulation  (infinitesimal  deformation  formulation)  with  the  following  Lame  constants:  A  = 
576.92  and  G  =  384.62. 

The  proposed  scheme  has  been  used  with  the  following  penalty  parameters:  kn  =  10^  and 
kt  =  10^  with  6  =  0.5. 
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Figure  8.15:  Forging  of  an  elastic  block  against  a  rigid  foundation.  Deformed  mesh. 

The  deformed  block  is  shown  in  Figure  8.15.  We  compare  the  results  of  our  scheme  with  the 
one  presented  in  Oden  &:  Pires[26],  for  the  case  with  fi  =  0.5  and  e  =  10~^,  where  €  = 

A  comparison  of  the  stresses  along  the  base  of  the  block  can  be  seen  in  Figure  8.16.  We  note 
that  there  is  good  agreement  between  the  curves  and  the  data  points. 

This  numerical  example  proves  not  only  that  the  scheme  has  good  stability  properties  but 
also  that  the  simulation  results  are  very  similar  to  the  results  of  other  implicit  schemes  which 
lack  these  properties. 

8.3.3  Oblique  impact  of  two  infinite  blocks 

We  present  in  this  section  the  case  of  an  oblique  impact  in  plane  strain  between  two  linearly 
elastic  blocks,  that  is,  we  use  an  elastic  formulation  that  considers  infinitesimal  deformations. 
The  top  block  has  an  initial  velocity  vq  =  (-10,  -10)  and  the  other  has  its  bottom  edge 
fixed.  Figure  8.17  shows  a  schematic  drawing  of  the  problem.  At  t  =  0.0,  the  block  is  resting 
on  the  bottom  block  (in  the  drawing  we  separated  the  blocks  for  clarity) . 

The  penalty  parameters  used  are  as  follows:  =  kt  =  10'^.  The  problem  has  been 

simulated  with  the  following  frictional  conditions:  fx  =  0.0  (frictionless)  and  n  =  0.4.  The 
displacement  of  point  A  (see  Figure  8.17)  in  both  directions  is  compared  with  the  results 
obtained  by  Chen  &  Yeh[8]. 

Figure  8.18  displays  the  horizontal  and  vertical  displacements  of  point  A  vs.  time,  for  the 
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Figure  8.16:  Forging  of  an  elastic  block  against  a  rigid  foundation.  Stress  curves  along  the 
base  of  the  block. 


5 


Figure  8.17:  Oblique  impact  of  two  elastic  blocks.  Schematic  drawing  of  the  problem. 
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Figure  8.18:  Oblique  impact  of  two  elastic  blocks.  Displacement  of  point  A. 

frictionless  and  frictional  cases.  As  expected,  the  horizontal  displacements  are  significantly 
reduced  by  the  friction  phenomenon,  while  the  vertical  displacements  on  the  rebound  increase 
when  friction  is  present. 

The  results  obtained  by  the  present  scheme  compare  well  with  the  results  from  [8].  In  Figures 
8.19,  8.20  and  8.21  we  show  the  deformed  mesh,  Cxx  and  Oyy,  respectively,  at  time  t  =  0.1 
for  the  frictional  case. 
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Figure  8.19:  Oblique  impact  of  two  elastic  blocks.  Deformed  mesh  at  t  =  0.1. 
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Figure  8.20:  Oblique  impact  of  two  elastic  blocks.  Stress  at  t  =  0.1  for  the  frictional 
case. 
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Figure  8.21:  Oblique  impact  of  two  elastic  blocks.  Stress  ayy  at  t  =  0.1  for  the  frictional 
case. 

8.4  Multibody  contact 

8.4.1  Impact  of  9  elcistic  disks 

We  show  the  numerical  results  of  the  impact  of  9  disks  enclosed  within  four  walls  (see  Figure 
8.22.  The  disks  have  Lame  constants  A  =  130.0  and  G  =  43.33,  and  density  p  =  8.93.  The 
disk  in  the  lower  left  corner  is  given  an  initial  velocity  vo  =  (0.5, 0.5)  and  the  the  disk  in  the 
upper  right  corner  is  given  an  initial  velocity  vq  =  (-0.5,  -0.5).  The  evolution  of  the  energy 
of  the  system  is  shown  in  Figure  8.23. 

During  the  simulation,  the  energy  calculated  using  the  mid-point  rule  increases  rapidly  be¬ 
cause  of  the  many  collisions.  In  contrast,  the  energy  calculated  using  the  proposed  scheme 
is  restored  to  its  original  value  after  each  collision.  Also,  the  simulation  using  the  mid-point 
rule  blows  up. 

The  motion  of  the  bodies  may  be  observed  in  Figure  8.22.  The  choice  of  material  properties 
and  initial  velocity  of  the  system  causes  the  bodies  to  suffer  significant  deformations.  There¬ 
fore,  this  example  shows  the  implementation  using  the  proposed  scheme  to  be  effective  in 
the  simulation  of  multiple  collisions  among  bodies  undergoing  finite  deformations  and  large 
rotations. 
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8.4.2  Impact  of  49  quasi-rigid  disks 

In  this  numerical  example  we  show  simulation  results  of  the  impact  of  49  elastic  disks  enclosed 
within  four  rigid  walls.  The  disks  have  Lame  constants  A  =  2000.0  and  G  =  1000.0,  and 
density  p  =  1.0.  The  left-most  column  of  disks  is  given  an  initial  velocity  vq  =  (0.5, 0.1). 
Figures  8.4.2  and  8.4.2  show  the  evolution  of  the  system.  Notice  that  clusters  of  bodies  are 
formed  at  various  instants,  i.e.  a  particular  body  may  be  in  contact  with  many  others  at 
certain  time  steps. 
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t  =  30.  t  =  35.  t  =  40. 


I 

I 


Figure  8.24:  Impact  of  49  quasi-rigid  disks.  Evolution  of  the  system, 


t  =  45. 


t  =  60. 


t  =  75. 


Figure  8.25:  Impact  of  49 


8.4.  MULTIBODY  CONTACT 


t  =  50. 


t  =  55. 


t  =  65.  t  =  70. 


t  =  80.  t  =  85. 

juasi-rigid  disks.  Evolution  of  the  system  (continued). 
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Chapter  9 


Conclusions 


9.1  Closure 

In  this  dissertation,  the  formulation  of  a  new  implicit  stable  time-stepping  algorithm  for 
dynamic  contact  problems  has  been  presented.  The  salient  points  of  the  proposed  scheme 
are  the  following: 

1.  In  the  frictionless  case,  the  new  contact  scheme  unconditionally  inherits  the  conser¬ 
vation  properties  of  the  underlying  continuum  problem,  thus  simulating  the  correct 
physical  system.  We  have  also  shown,  by  means  of  numerical  examples,  that  the  re¬ 
sults  obtained  may  be  significantly  different  if  non-conserving  schemes  are  used,  even 
if  no  numerical  blow-up  occurs.  With  almost  no  modifications,  our  scheme  can  also  be 
used  also  in  quasi-static  simulations. 

2.  In  the  frictional  case,  the  scheme  is  proven  to  be  unconditionally  dissipative  thus  imi¬ 
tating  the  real  physical  system.  It  has  been  shown,  by  means  of  numerical  examples, 
that  standard  frictional  schemes  may  predict  energy  increases,  thus  yielding  unphys¬ 
ical  behaviour.  Furthermore,  our  frictional  scheme  can  also  be  used  in  quasi-static 
problems. 
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3.  Using  penalty  methods,  we  are  able  to  enforce  the  constraint  on  the  gap  velocity,  in  ad¬ 
dition  to  satisfying  efficiently  the  usual  impenetrability  constraint.  This  enhancement 
does  not  affect  the  previously  mentioned  conserving  properties  of  the  proposed  scheme. 
The  enforcement  of  this  constraint  enables  the  proposed  scheme  to  circumvent  some 
of  the  difficulties  present  in  short-duration  simulations  in  which  oscillatory  behaviour 
tends  to  overwhelm  the  real  solution. 

4.  Since  the  energy  is  under  control  at  all  times,  we  have  an  (energy)  stable  algorithm. 
This  property  yields  a  robust  scheme  which  is  easily  implementable,  with  only  slight 
changes,  into  the  standard  penalty  formulations. 

5.  We  have  also  discussed  the  introduction  of  positive  high  frequency  energy  dissipation 
for  short-duration  simulations.  The  amount  of  dissipation  can  be  regulated  depending 
on  the  problem  at  hand,  so  as  to  eliminate  the  desired  degree  of  oscillation  in  the 
results.  This  enhancement  requires  only  a  slight  modification  of  the  original  energy 
restoring  scheme. 

6.  The  implementation  of  the  proposed  scheme  in  a  multibody  formulation  has  been 
described  in  detail.  Several  concepts  from  the  discrete  element  method  were  found 
to  be  versatile  and  adaptable  to  a  finite  element  contact  formulation,  in  particular 
the  contact  detection  algorithms.  The  introduction  of  sorting  procedures  prior  to  the 
contact  resolution  phase,  namely  the  closest  point  procedure,  proved  to  be  eflicient 
when  dealing  with  multiple  bodies. 

7.  Several  numerical  examples  have  been  presented  to  show  improved  results  when  using 
the  new  scheme,  demonstrating  the  value  of  satisfying  the  conservation  laws  of  the 
dynamical  system  in  question.  The  computational  effort  to  simulate  these  problems 
was  not  higher  when  using  our  scheme  than  when  using  more  conventional  schemes. 

f  y.  ■ «  nr-f  ■  ■ 
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9.2  Future  work 

Contact  problems  are  in  fact  coupled  problems,  in  the  sense  that  there  are  multiple  bodies 
which  can  be  simulated  independently  unless  they  come  in  contact  with  one  another.  In  a 
dynamical  system,  keeping  track  of  the  bodies  which  are  in  contact  can  become  cumbersome. 
Therefore,  the  use  of  staggered  algorithms  may  work  towards  improving  the  efficiency  of  the 
existing  algorithms.  Of  course,  it  is  desirable  that  these  staggered  schemes  have  the  same 
stability  properties  as  the  coupled  scheme. 

Prom  a  practical  standpoint,  adding  other  physical  phenomena  to  the  contact  problem  is 
an  attractive  option.  Some  examples  are  the  simulation  of  heat  during  contact,  where  the 
source  of  heat  may  be  friction  or  plastic  deformation. 

Future  work  may  also  concentrate  on  the  development  of  more  efficient  contact  detection  al¬ 
gorithms  which  make  use  of  object-oriented  programming,  especially  in  the  area  of  multibody 
simulations. 
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Appendix  A 


Consistent  Linearization  of  the 
Proposed  Schemes 


For  two  dimensional  problems,  we  develop  in  this  appendix  the  consistent  linearization  of 
the  time  stepping  schemes  we  have  presented  in  previous  chapters. 

In  Section  A.l,  the  linearized  equations  of  the  problem  are  developed  and,  in  Section  A.2 
and  A.3  the  linearized  contributions  of  the  normal  and  frictional  contact  forces  are  presented 
respectively. 


A.l  The  linearized  equations 

Consider  the  system  of  ordinary  differential  equations  4.42  and  4.43  and  introduce  in  them 
the  definition  of  the  nodal  momenta  p  defined  in  equation  4.44.  In  this  way  the  linearization 
of  equation  4.43  remains  unchanged  even  with  the  introduction  of  the  modified  nodal  linear 
momenta  by  means  of  the  mass  penalty  rrip  for  the  nodes  in  contact.  Equation  4.43  becomes 
nonlinear  becomes  nonlinear  due  to  this  penalization  thus  we  need  to  consider  its  linearization 
too.  Therefore,  we  consider  the  following  residuals: 
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Using  the  Newton-Raphson  methodology,  we  perform  the  consistent  linearization  of  the 
residuals  A.I  and  A.2,  yielding 
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are  the  updates  between  iterations  (i)  and  (i  +  1)  in  time  step  [tn,<n+i]- 


Consider  the  following  notation 


(A.7) 


with  =  iAd^'^i^  for  the  contributions  of  the  internal  elastic  forces  to  the  tangent 

stiffness,  and 


(A.8) 


for  the  contribution  of  the  contact  forces.  . 
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Using  equation  A.l,  we  eliminate  from  equation  A.2.  The  introduction  of  equations 

A. 7  and  A.8,  leads  to 


-  -H 


(i) 

'd 


(A.9) 


The  previous  equation  yields  Ad^+P,  which  in  turn  is  used  to  calculate  Ap^+p.  We  are  able 
to  find  the  expression  for  the  update  of  the  nodal  velocities  using  the  definition  of  the 
modified  nodal  momenta  given  in  equation  5.38  and  is  given  by 
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S=1 


=:  M*(*+0 


Notice  that  the  update  equation  A.  10  is  only  applied  to  the  nodes  participating  in  the 
mass  penalty  regularization.  For  the  rest  of  the  nodes,  equation  A.l  is  linear  so  that  its 
linearization  leads  to 


=  ’’i’ii  +  ■  (A.ll) 

Remark  A.l  An  implementation  avoiding  the  use  of  nodal  momenta  p  for  the  nodes  in 
contact  can  be  easily  devised  by  considering  the  linearized  version  of  equation  A.  10.  Details 
are  omitted. 


A. 2  The  contact  stiffness 

The  linearization  of  the  contact  force  fc,  defined  by  equation  4.36  as 


/\  JS,C  J 
5=1 


(A.12) 
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with 
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is  given  by: 


with 


A/. 


(n+i) 


2'  _ 


'^slave 

A 

5=1 


(A.13) 


A/1ri'  = 
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material{normal  part) 


^  AtrHg,„+i^  + 

5eomefrzc(norma/  part)  material{tangential  part)  geometric(tangential  part) 


(A.14) 


and  we  have  defined  the  expression  Gs,t  in  equation  4.38  and  the  expression  Hg^t  in  6.4.  We 
use  the  following  notation  to  rewrite  equation  A.14: 


A^^^^  =  KcAdn+i  (A.15) 

=  +  Kf^°  +  irf  +  fcf "  .  (A.16) 

Next  we  develop  the  linearizations  for  the  normal  and  tangential  contributions  of  the  contact 
force  to  the  overall  consistent  stiffness  tangent  matrix,  that  is  K,. 

A. 2.1  The  normal  contact  stiffness 

The  material  part  of  the  normal  cntact  stiffness  is  obtained  through  the  linearization  of  the 
normal  contact  pressure  and  is  given  by 
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if  4n+l  ^  din  > 
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Atyv  =  < 


(A.17) 


A.2.  THE  CONTACT  STIFFNESS 


(A.22) 


The  expression  refers  to  the  real  gap  calculated  through  the  closest  point  projection 

at  the  configuration  at  is  the  tangent  vector  to  the  master  surface  at  the  point 

of  contact  ( i.e.  =  0)  and  Ig  is  the  length  of  the  surface  element  of  the  master 

surface  corresponding  to  the  given  slave  node  S. 

Let  us  define  the  following  auxiliary  vectors; 
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The  geometrical  part  of  the  tangent,  arising  from  the  change  of  normal  and  contact  point  in 
is  as  follows 
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after  an  involved  calculation  found  in  Laursen  [18].  The  final  expressions  for  the  contact 
stiffness  are  then  given  by 


-r^TflCbt 

iljv 


U'  {gin+l)  - 


ffs,n+l 


9s, n 


^5,n+i 


^s,n+i 


Cl  ^  C2 

_n  1 - Z'T'  , 

2  -^5,72+^  2 


(A.26) 


and 
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with  the  difference  quotient  in  the  first  term  of  replaced  by  U"  if  9s,n+i  =  9s, n- 

We  note  the  non-symmetry  of  the  material  part  as  it  occurs  with  its  counterpart  for  the 
energy-momentum  conserving  algorithms  considered  in  this  paper  for  the  continuum. 


A. 2. 2  The  tangent  contact  stiffness 

The  material  part  of  the  tangent  contact  stiffness  is  obtained  through  the  linearization  of 
the  tangent  friction  traction  ty.  The  expression  for  tx  depends  on  whether  we  have  a  stick 
or  a  slip  condition. 
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1.  Stick  phase 


Atx  - 

=  KrMu^Ae^+i  ,  (A.28) 

where 

=  »=,.+iA5.+i  +  AhJ„^i(X.+,  - 2.)  ,  (A.29) 

and  we  can  recognize  that 
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Using  equation  A.29  and  A.30,  we  can  rewrite  the  former  as 
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Finally,  using  the  previous  algebraic  manipulations,  we  obtain 
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2.  Slip  phase 
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where  we  know  from  the  previous  section  that 
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with  the  definition  the  parameters  Ci  and  C2,  and  the  arrays  Ts,t  and  Ds,t,  given  above. 
The  expression  in  parentheses  in  equation  A.34  is  replaced  by  if  g^^n+i  =  9s  n- 

Also, 


ftrial 

11  ^ 


where  we  denote  by  Ss,t  the  following  array: 
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Thus,  we  can  express  the  material  contribution  of  the  friction  force  to  the  overall 
consistent  tangent  matrix  as  follows: 
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For  the  geometric  part  of  the  matrix,  we  linearize  the  expression  developed  in 

Laursen[18],  as  follows: 
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where  we  have  defined  as 
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Thus,  the  final  expression  for  the 
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tangent  matrix  is: 
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